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Research  is  needed  to  determine  the  microscopic  origins  of  material  instability.  A 
microscopic  investigation  becomes  the  rational  means  to  examine  shear  band  instability 
which  has  a  thickness  of  several  grains  (i.e.,  microscopic  dimension). 

Research  objectives  and  report  organization 

Many  continuum  theories  and  numerical  simulations  have  been  proposed  to  describe 
shear  band  instability;  however  they  rely  on  semi-intuitive  arguments  not  always 
founded  on  physical  observations  on  particulate  media.  The  main  research  objective  is  to 
investigate  the  physical  origins  of  shear  band  instability  in  particulate  media,  and  the 
assumptions  of  higher-order  continuum  theories  proposed  to  account  for  instabilities  in 
granular  media.  The  particular  research  objectives  were  (1)  to  review  the  work  in 
computational  granular  mechanics  relevant  to  the  formation  of  shear  bands  in  granular 
media,  (2)  to  explore  a  critical  micro-macro  transition  relevant  to  material  instability  (i.e., 
stress  symmetry),  and  (3)  to  devise  laboratory  experiments  for  generating  useful 
experimental  data  sets  relevant  to  material  instability. 

This  document  is  organized  in  three  sections.  Following  the  introduction  and  research 
background,  the  first  section  reviews  past  work  on  computational  granular  mechanics. 
The  second  section  investigates  the  transition  between  discrete  and  continuous  fields,  and 
especially  the  symmetry  of  stress  tensors  in  granular  media.  The  third  and  last  section 
summarizes  the  experimental  work  on  discrete  models. 
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PART  I.  GRANULAR  MECHANICS 


Introduction 

The  first  discrete  modeling  of  soils  can  be  traced  to  Hertz  (1882)  who  formulated  a 
contact  law  between  spheres,  and  Reynolds  (1885)  who  proposed  a  dilatancy  theory. 
Dantu  (1957)  and  Schneebli  (1955)  idealized  real  soils  as  assemblies  of  rigid  rods,  and 
noticed  some  striking  similarities  between  the  mechanical  responses  of  these  mechanical 
analogs  and  real  soils.  Duffy  and  Mindlin  (1957),  Deresiewicz  (1958),  and  Thurston  and 
Deresiewicz  (1959)  examined  the  response  of  soil  models  made  of  spheres.  Biarez  (1962) 
used  glass  beads  and  duralumium  rods  to  examine  the  elastic  and  limit  response  of  soils, 
and  applied  his  observations  to  analyze  practical  problems  in  geotechnical  engineering. 
These  pioneer  works  were  later  followed  by  photoelastic  investigations  (e.g.,  Drescher, 
1976;  and  Descher  and  Josselin  de  Jong,  1972)  to  visualize  stresses  within  granular 
media. 

The  discrete  modeling  of  soils  benefited  substantially  from  the  development  of  computers 
in  the  1970's.  The  computational  discrete  modeling  of  soils  can  be  attributed  to  Cundall, 
who  developed  the  computer  code  BALL  (Cundall  and  Strack,  1978-79).  At  this  time, 
computers  had  very  limited  capabilities,  comparable  to  those  of  today's  hand-held 
calculators.  They  were  slow,  and  had  limited  memory  and  storage  capacity.  Yet,  Cundall 
developed  BALL,  a  program  which  many  researchers  still  use.  The  program  is  fully 
documented  in  a  two-volume  report  to  the  US  National  Science  Foundation  (Cundall  and 
Strack,  1978-79).  Since  1978,  many  researchers  have  adapted  the  original  version  of 
BALL  to  solve  specific  problems.  Researchers  have  realized  the  power  of  computer 
simulations  to  understand  the  mechanics  of  granular  materials.  The  major  advantage  of 
computations  over  real  experiments  is  the  generation  of  abundant  information  on  particle 
displacements,  contact  forces  and  other  physical  quantities,  which  can  be  processed 
rapidly  to  comprehend  the  physics  of  granular  assemblages. 

This  chapter  reviews  the  basics  of  the  discrete  modeling  of  granular  media,  with  the 
intent  of  providing  readers  with  some  understanding  of  granular  material  behavior  based 


on  the  first  principles  in  mechanics  and  computational  methods.  Three  particular  aspects 
of  computational  granular  mechanics  will  be  covered:  (1)  physical  modeling  of  granular 
media,  (2)  numerical  methods  for  computational  granular  mechanics,  and  (3)  transitions 
from  discrete  to  continuous  media.  Before  going  into  these  subjects,  we  will  first  try  to 
answer  the  following  questions.  What  exactly  is  discrete  modeling?  What  are  the 
applications  of  discrete  modeling?  Has  discrete  modeling  any  advantages  over  continuum 
mechanics?  What  are  the  limitations  of  discrete  modeling?  What  are  the  main 
computational  tools  in  discrete  modeling? 

Examples  of  discrete  modeling  in  soil  mechanics 

Discrete  modeling  can  best  be  described  by  considering  two  particular  examples:  one  is 
relevant  to  the  study  of  constitutive  behavior,  and  the  other  to  the  failure  of  foundations. 

Example  1 

As  shown  in  Fig.  1,  the  sample  was  constructed  with  1848  cylindrical  rods  of  diameter  4, 
6  and  8  mm.  There  are  1040  particles  of  4-mm  diameter,  532  of  6-mm  diameter,  and  266 
of  8-mm  diameter.  The  sample  slenderness  ratio  is  2.43.  The  rods  were  made  of  identical 
transparent  acrylic  material,  which  has  an  average  density  of  1.30  g/cm3,  and  were  cut  to 
a  10-cm  length.  Their  front  ends  were  half  painted  in  black  to  visualize  their  rotation  and 
displacement  simultaneously.  Figures  1  and  2  show  the  experimental  setup.  The  sample  is 
enclosed  in  a  0.1 5 -mm  thick  transparent  latex  membrane  with  two  circular  loading 
platens  at  the  top  and  bottom.  The  membrane  is  clamped  to  the  platens  with  compression 
rings.  During  the  sample  fabrication,  the  membrane  is  stretched  on  a  rectangular  mold, 
and  the  rods  are  manually  positioned  inside  the  mold  with  their  axis  parallel  to  one 
another.  During  the  test,  the  specimen  is  subjected  to  a  vacuum  inside  the  membrane, 
which  is  equivalent  to  applying  an  external  constant  confining  pressure  equal  to  95  kPa. 
The  axial  compression  is  applied  by  raising  the  lower  platen  at  the  constant  rate  of  5.6 
mm/min,  while  the  upper  platen  remains  fixed.  Figure  2  shows  a  side  view  of  the 
experimental  setup.  The  photograph  of  Fig.  1  was  obtained  by  using  a  35-mm  camera 
positioned  approximately  one  meter  in  front  of  the  sample,  with  its  aim  perpendicular  to 
the  front  of  the  sample.  The  particles  were  lighted  from  behind  by  using  a  light  box.  The 


-  Page  7  - 


sample  was  tested  by  increasing  the  axial  strain  while  keeping  the  confining  pressure 
constant.  The  stress-strain  curve  is  shown  in  Fig.  3.  The  axial  strain  e  is 


£=  Ah/ho  (1*0 

where  Ah  is  the  vertical  displacement  of  the  lower  platen,  and  h0  is  the  initial  sample 
height.  The  axial  stress  a  is 

a=F/A0  O-2) 

where  F  is  the  measured  axial  load,  and  A0  is  the  initial  cross-sectional  area.  The  Initial 
Young's  modulus  was  300  MPa,  and  the  peak  friction  angle  was  18.4  deg. 
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Figure  3.  Stress-strain  response  of  sample  shown  in  Fig.  1  (the  confining  pressure  is  95 
kPa). 

Example  2 

Figure  4  shows  an  early  example  of  the  application  of  discrete  modeling  for 
understanding  the  failure  mechanism  under  shallow  footing  (Lambe  and  Whitman,  1979). 
This  photograph  was  taken  by  using  a  slow  shutter  speed,  which  blurs  the  moving 
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particles.  It  shows  the  pattern  of  motion  at  failure  within  a  stack  of  duralumium  rods 
loaded  by  a  rigid  punch.  The  rods  are  150-mm  long  and  are  of  two  shapes  and  sizes 
(round,  3  mm  and  6  mm  diameter;  and  hexagonal  4.8  and  7,9  mm  across  flats)  to 
simulate  the  interlocking  which  occurs  in  actual  soils.  Additional  examples  of  discrete 
models  for  understanding  the  failures  of  foundations  and  retaining  walls  can  be  found  in 
(Biarez,  1962). 


Figure  4.  Failure  zone  under  a  shallow  foundation  (Lambe  and  Whitman,  1979). 

Applications  of  discrete  modeling  in  engineering  and  applied  sciences 

One  may  attempt  to  survey  the  applications  of  discrete  modeling  in  engineering  and 
applied  sciences  by  consulting  the  proceedings  of  the  following  conferences: 

(1)  the  1st  International  Conference  on  Powders  and  Grains  (Biarez  and  Gourves,  1989); 

(2)  the  2nd  International  Conference  on  Powders  and  Grains  (Thornton,  1993); 

(3)  the  3rd  International  Conference  on  Powders  and  Grains  (Behringer  and  Jenkins, 
1997); 

(4)  the  1st  US  Conference  on  Discrete  Element  Methods  (Mustoe  et  al.,  1989);  and 

(5)  the  2nd  International  Conference  on  Powders  and  Grains  (Wiliams  and  Mustoe, 
1993). 
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Discrete  modeling  is  used  for  many  applications  in  engineering  and  applied  sciences.  A 
review  of  the  literature  on  discrete  modeling,  which  is  by  no  means  exhaustive,  indicates 
that  past  work  can  be  sorted  in  the  following  general  categories: 

•  Experimental  studies  of  granular  media :  (e.g.,  Biarez,  1962;  Dantu,  1957;  Drescher, 
1976;  Drescher  and  Josselin  de  Jong,  1972;  Meakin  and  Skeltorp,  1993;  Nakase  et  al., 
1992;  Oda,  1972-1973;  Subbash  et  al.,  1991;  Sukla  and  Rossmanith,  1982;  and  Supel, 
1985). 

•  Studies  of  contact  between  particles:  (e.g.,  Anandarajah,  1994;  Azarkhin,  1988; 
Bardet  and  Huang,  1992;  Bentall  and  Johnson,  1967;  Carter,  1926;  Deresiewicz, 

1958;  Domenech  et  al.,  1987;  Fabrikant,  1986, 1988;  Foerster  et  al.,  1994;  Goodman, 
1962;  Hertz,  1882;  Jagota  et  a.,  1997;  Johnson  et  al.,  1997;  Johnson,  1958, 1985; 
Kalker,  1979,  Lecomu,  1905;  Lee,  1966;  Lesbur  et  al.,  1997;  Mindlin,  1949,  Mindlin 
and  Deresiewicz,  1953;  Misra,  1995;  Reynolds,  1885, 1895;  tabor,  1955;  Vermeulen 
and  Johnson,  1964;  Wells,  1997;  Witters  and  Duymelinck,  1986). 

•  Numerical  techniques  for  discrete  modeling:  (e.g.,  Bardet  and  Proubet,  1991;  Bashir 
and  Goddard,  1991;  Borja  and  Wren,  1995;  Chang  and  Misra,  1989;  Cundall,  1988, 
1989;  Cundall  and  Strack,  1979;  Ghaboussi  and  Barbosa,  1990;  Goddard  et  al.,  1993; 
Hahn,  1988;  Hart  et  al.,  1988;  Jean,  1995;  Jean  and  Moreau,  1996;  Kishino,  1988; 
LaBudde  and  Greenspan,  1976;  Meakin  and  Skeltorp,  1993;  Moreau,  1966, 1988, 
1994,  1995;  Mujinza  et  al,  1993;  O'connor  et  al.,  1993;  Papadrakakis,  1981;  Park  and 
Underwood,  1980;  Ting,  1992;  Ting  and  Corkum,  1992;  Ting  et  al.,  1989,  1993; 

Trent  and  Margolin,  1992;  Underwood,  1983;  Underwood  and  Park,  1980;  Walton, 
1980, 1983;  Walton  and  Braun,  1986;  Walton  et  al.,  1988;  Williams  and  Mustoe, 
1987;  Zhuang  and  Goddard,  1993;  Zhuang  et  al.,  1995). 

•  Elastic  behavior  of  granular  media:  (e.g.,  Bathurst  and  Rothenburg,  1988;  Chang  et 
al.,  1990;  Changet  al.,  1995;  Dubujet  etal.,  1997;  Home,  1965,  1969;  Johnson  et  al., 
1997). 
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Failure  of  granular  media :  (e.g.,  Chang,  1993;  Sun  and  Thornton,  1994;  Thornton, 
1979). 


•  Stress-strain  behavior  of  granular  media:  (e.g.,  Bardet,  1994;  Bathurst  and 
Rothenburg,  1990;  Chang,  1993;  Chang  and  Liao,  1990;  Chang  and  Misra,  1990; 
Chang  et  al.,  1992;  Christoffersen  et  al.,  1981;  Cundall  et  al.,  1982;  Daudon  et  al., 
1997;  Deresiewicz,  1958;  Dobry  and  Ng,  1992;  Duffy  and  Mindlin,  1957;  Houlsby, 
1981;  Ke  and  Bray,  1995;  Koenders,  1987;  Krawietz  1982;  Laalai  et  al.,  1995;  Lun 
and  Bent,  1993;  Matsuoka  and  Yamamoto,  1994;  Mehrabadi  et  al.,  1993;  Nemat- 
Nasser  and  Mehrabadi,  1984;  Ng  and  Dobry,  1994;  Ng,  1992;  Reynolds,  1895; 
Rothenburg  and  Bathurst,  1989, 1991, 1992;  Rowe,  1962;  Tatsuoka  et  al.,  1990; 
Thornton,  1979, 1994;  Thornton  and  Randall,  1988;  Thornton  and  Sun,  1994; 
Thurston  and  Deresiewicz,  1959;  Weber,  1996). 

•  Higher-order  continuum  theories: 

•  Second  order  gradient  continuum:  (e.g.,  Chang  and  Gao,  1995). 

•  Micropolar  continuum:  (e.g.,  Bardet  and  Proubet,  1992;  Bardet  and  Huang,  1992; 
Chang  and  Ma,  1991;  Diepolder  et  al.,  1991;  Kanatani,  1979;  Sternberg,  1968). 

•  Strain  localization  and  shear  bands:  (e.g.,  Bardet  and  Proubet,  1991,  1992;  Cundall, 
1989,  Desrues,  1984,  Desrues  and  Duthilleul,  1984;  Desrues  et  al.,  1985;  Moreau, 
1996;  Muhlhaus  and  Vardoulakis,  1987;  Nakase  et  al.,  1992;  Oda,  1993;  Oda  et  al., 
1997). 

•  Powders  and  sintering  processes:  (e.g.,  Aizawa  et  al.,  1993;  Greening  et  al.,  1997; 
Hong,  1997;  Lian  et  al.,  1997;  Tamura  and  Aizawa,  1993). 

•  Suspension:  (e.g.,  Goddard,  1977,  1986). 

•  Granular  Flow:  (e.g.,  Foerster  et  al.,  1994;  Campbell  and  Brennen,  1985;  Savage  and 
Jeffrey,  1981). 
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•  Fluid  and  solid  mixtures,  including  fluidized  beds :  (e.g.,  Kawaguchi  et  al.,  1995; 

Sudji  et  al,  1992, 1993;  Tsuji,  1997). 

•  Rock  mechanics’,  (e.g.,  Bardetand  Scott,  1985;  Cundall,  1971, 1981,  1988;  Shi,  1993; 
Shi  and  Goodman,  1988, 1989). 

•  Physics  of  granular  media:  (e.g.,  Clement  et  al.,  1992;  Knight  et  al.,  1993;  Manna 
and  Herrmann,  1991;  Meakin  and  Skeltorp,  1993;  Meftah  et  al.,  1993;  Radjai  et  al., 
1996). 

Discrete  and  continuous  modeling 

Continuum  mechanics  is  a  powerful  approach  to  solve  scientific  and  engineering 
problems.  It  is  based  on  mathematical  assumptions,  which  are  well  described  in  the 
continuum  mechanics  literature  (e.g.,  Eringen,  1967;  Truesdell,  1985).  It  formulates 
engineering  problems  as  mathematical  boundary  value  problems  (BVP).  The  main 
components  are  the  governing  equilibrium  equations  (i.e.,  partial  differential  equations 
translating  basic  physical  balances,  such  as  stress  equilibrium);  the  boundary  conditions 
(prescribed  values  of  the  unknown  quantities  or  their  derivatives  on  external  surfaces); 
and  the  constitutive  equations  (generalized  relation  between  stress  and  strain,  or  their 
respective  rates).  In  the  past  20  years,  the  continuum  approach  has  been  implemented 
using  finite  elements  and  finite  differences,  and  has  successfully  been  applied  to  solve 
engineering  problems  (e.g.,  Zienkiewicz  and  Taylor,  1991). 

One  of  the  major  assumptions  of  continuum  mechanics  is  that  the  material  properties  can 
be  scaled  from  small  laboratory  samples  to  large  material  masses  using  constitutive 
relations.  Researchers  have  now  produced  a  myriad  of  constitutive  relations  for  various 
types  of  engineering  materials.  Many  constitutive  models  are  clever  fittings  of 
experimental  results  in  the  laboratory.  Unfortunately,  most  models  are  not  based  on 
physics. 

Through  the  use  of  bifurcation  theory,  the  continuum  approach  was  discovered  to  exhibit 
problems,  especially  associated  with  the  loss  of  uniqueness  of  boundary  value  problems. 
The  problems  mainly  result  from  local  material  instability  and  limitations  of  constitutive 
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equations  (e.g.,  Vardoulakis  and  Graf,  1985).  These  findings  questioned  some 
applications  of  continuum  mechanics  to  soil  mechanics,  and  pointed  out  the  need  for 
seeking  a  fundamental  understanding  of  material  behavior.  Discrete  modeling  can  be  of 
great  help  to  continuum  mechanics,  for  not  only  developing  constitutive  models  based  on 
physics,  but  also  for  understanding  the  physical  origins  of  material  instability,  and  the 
limitations  of  continuum  mechanics. 

Limitations  of  discrete  modeling  in  soil  mechanics 

Discrete  modeling  has  obvious  limitations  in  soil  mechanics.  The  sheer  number  of 
individual  soil  particles,  especially  those  with  smaller  diameters,  within  soil  masses  of 
practical  interest  to  engineering,  prohibits  the  simulation  of  their  overall  response  even 
with  the  most  advanced  computers.  For  instance,  in  a  cubic  centimeter  of  soil,  there  may 
be  as  many  as  5xl012  clay  particles,  when  those  are  assumed  to  be  identical  square 
platelets  (1-pm  x  1-pm  x  0.1  -|im),  and  the  void  ratio  is  equal  to  one.  Therefore,  the 
largest  computers  yet  built  would  not  be  sufficient  to  handle  1  cm3  of  fine-grained  soils, 
which  is  quite  irrelevant  for  engineering  purposes.  The  number  of  soil  particles  in  a 
volume  decreases  roughly  with  the  cube  of  the  particle  size.  In  a  cubic  centimeter,  there 
may  be  as  many  as  1,000  particles  of  coarse  sand,  when  those  are  assumed  rounded  with 
1-mm  diameter  and  the  void  ratio  is  still  equal  to  one.  Yet  again,  the  most  powerful 
computers  of  today,  with  numerous  parallel  processors,  would  have  extreme  difficulty 
handling  1  m3  of  coarse  sand,  which  corresponds  to  lxlO9  particles. 

In  view  of  the  excessively  large  number  of  particles  in  soils,  it  seems  more  feasible  to 
combine  the  advantages  of  discrete  modeling  and  continuum  mechanics  to  solve 
engineering  problems.  In  summary,  discrete  and  continuum  mechanics  should  be 
perceived  as  complementary,  not  adversary,  tools  in  soil  mechanics  to  understand  the 
mechanics  and  physics  of  soil  behavior. 

Numerical  methods  for  discrete  modeling 

The  main  numerical  tool  of  discrete  modeling  is  the  discrete  element  method  (Cundall 
and  Strack,  1979).  Cundall  (1989)  proposed  that  this  appellation  apply  to  a  computer 
program  only  if  it  (a)  allows  finite  displacements  and  rotations  of  discrete  bodies, 
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including  detachment,  and  (b)  recognizes  new  contacts  automatically  as  the  calculation 
progresses.  To  my  knowledge,  there  are  eight  main  classes  of  numerical  methods 
corresponding  to  this  definition. 

1 .  Distinct  element  methods  (DEM).  These  methods  use  explicit  and  time-marching 
algorithms  to  solve  the  equation  of  motion.  Bodies  may  be  rigid  or  deformable,  but 
contacts  are  always  deformable.  Examples  of  DEM  codes  are  TRUBAL  (Cundall, 
1988);  UDEC  (Cundall,  1980);  3DEC  (Hart  et  al.,  1988);  DIBS  (Walton,  1980); 
3DSHEAR  (Walton  et  al.,  1988);  and  JP2  (Bardet  and  Proubet,  1989).  In  granular 
statics,  DEM  calculates  the  equilibrium  states  of  particle  systems  by  using  dynamic 
transitions,  the  convergence  of  which  are  generally  accelerated  and  optimized  by 
introducing  an  artificial  viscous  damping.  Such  optimizations  include  density  scaling 
(Cundalll,  1982)  and  adaptative  dynamic  relaxation  (Bardet  and  Proubet,  1991). 

2.  Modal  methods.  These  methods,  which  are  similar  to  DEM  in  the  case  of  rigid  bodies, 
use  modal  superposition  for  deformable  bodies,  (e.g.,  Williams  and  Mustoe,  1987).  A 
representative  code  is  CICE  (Hocking  et  al.,  1985). 

3.  Discontinuous  deformation  analysis  (DDA).  In  DDA  (Shi,  1993),  contacts  are  rigid, 
and  bodies  may  be  rigid  or  deformable.  The  condition  of  no-interpenetration  is 
achieved  by  an  iteration  scheme;  the  body  deformability  comes  from  superposition  of 
strain  modes. 

4.  Momentum-exchange  methods.  The  contacts  and  bodies  are  both  rigid:  momentum  is 
exchanged  between  two  contacting  bodies  during  an  instantaneous  collision. 

Frictional  sliding  can  be  represented  (e.g.,  Campbell  and  Brennen,  1985;  Hahn, 

1988). 

5.  Multibody  Dynamics  methods  (MDM).  Moreau  (1966)  treats  the  problem  of  non¬ 
penetrability  by  using  Convex  Analysis.  The  unilateral  mechanical  constraints  of 
frictional  non-penetrability  are  mathematically  formulated  in  (Brogliato,  1996; 
Delassus,  1917;  Moreau,  1966, 1988, 1993, 1994, 1995;  Pfeiffer  and  Glocker,  1966). 
MDM  was  implemented  in  (Jean,  1995)  using  an  implicit  algorithm,  and  was  applied 
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by  several  investigators  to  examine  various  aspects  of  particulate  behavior  (e.g., 
Daudon  et  al.,  1997).  MDM  considers  only  purely  rigid  bodies,  and  ignores  the 
deformability  of  individual  grains  and  contacts. 

6.  Structural  Mechanics  methods  (SMM).  These  methods  derive  from  the  numerical 
techniques  used  in  computational  mechanics,  especially  finite  element  methods  for 
continuum  plasticity  and  contact  mechanics  (e.g.,  Borja  and  Wren,  1995;  Kishino, 
1988;  Zienkiewicz  and  Taylor,  1991).  The  equilibrium  equations  for  the  system  of 
particles  are  solved  quasi-statically,  and  not  dynamically,  which  eliminates  the 
spurious  oscillations  of  dynamic  relaxation.  SMM  obtain  the  incremental  transition 
between  equilibrium  states  by  relying  on  a  tangential  stiffness  matrix,  which  is 
determined  from  the  local  contact  stiffness  between  particles.  Unfortunately,  this 
tangential  operator  consumes  a  large  amount  of  computer  memory,  rendering  it 
inapplicable  to  large  numbers  of  particles,  and  often  becomes  singular  and  causes 
numerical  problems.  However,  SMM  guarantees  a  strict  convergence,  when  the 
physical  problem  has  a  solution,  and  is  capable  of  detecting  bifurcation  points.  This 
approach,  which  benefits  from  the  progress  in  finite  element  methods  for  structural 
and  continuum  mechanics  (e.g.,  Bathe,  1996;  Zienkiewicz  and  Taylor,  1991),  reveals 
that  there  are  many  similarities  between  the  numerical  techniques  in  finite  and 
discrete  element  methods. 

7.  Mean  Field  method  (e.g.,  Bashir  and  Goddard,  1991;  Zhuang  et  al.,  1 995).  The 
transitions  between  static  states  are  calculated  by  imposing  a  mean  field  of 
displacement  and  rotation  to  the  particles,  and  by  restoring  a  new  equilibrium 
configuration  by  means  of  incremental  motions  or  "fluctuations"  of  each  particle 
about  the  mean.  The  fluctuating  motions  of  individual  particles  are  determined 
statically  by  a  global  stiffness  matrix,  which  is  determined  from  the  local  contact 
stiffness  between  particles,  as  in  structural  mechanics.  The  problems  arising  from  the 
matrix  singularity  are  overcome  by  solving  the  linear  equations  with  the  well-known 
relaxation  method  (Southwell,  1940). 


8.  Energy  minimization  method  (e.g.,  Chichili  et  al.,  1993).  The  transition  between  static 
states  is  obtained  by  incorporating  the  geometric  constraints  of  non-penetrability, 
Mohr-Coulomb  friction  criterion  for  relative  sliding,  and  a  minimization  function  that 
translates  the  energy  dissipated  at  the  contacts  due  to  internal  friction.  The 
interparticle  forces  are  calculated  at  each  time  step  by  using  an  explicit  finite 
difference  scheme  based  on  linear  programming  technique. 

Following  the  classification  in  (Cundall,  1989),  Fig.  5  summarizes  the  attributes, 
advantages  and  shortcomings  of  the  methods  listed  above,  which  include: 

•  contact  and  body  stiffness 

•  number  and  shapes  of  bodies 

•  capabilities  of  fracturing  individual  particles 

•  packing  density 

•  amplitude  of  displacement  and  strain,  and 

•  static  and  dynamic  capabilities. 

A  shown  in  Fig.  5,  the  method  performances  are  grouped  in  three  categories,  ranging 
from  good  to  not  applicable.  At  the  present,  it  is  difficult  to  conclude  on  the  superiority  of 
a  particular  method,  due  to  a  lack  of  in-depth  comparative  studies.  In  this  chapter,  we  will 
not  cover  all  these  methods,  but  will  only  introduce  the  basics  of  discrete  element 
methods. 

Physical  modeling  of  granular  media 

Geometry  of  grains 

As  described  in  soil  mechanics  textbooks  (e.g.,  Bardet,  1997;  Lambe  and  Whitman, 

1979),  soil  grains  have  irregular  and  various  shapes  including  spheres,  ellipsoids, 
platelets,  cylinders,  and  tubes,  when  they  are  observed  with  the  naked  eye  and 
microscopes.  Their  wide  range  of  grain  sizes  vary  from  colloids  (<1  pm)  to  boulders 
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(>100  cm).  The  diversity  of  grain  shape,  size,  distribution,  and  structures  are  one  of  the 
major  causes  of  the  multiplicity  of  soil  behavior  observed  in  the  laboratory  and  the  field. 

For  the  sake  of  simplification,  hereafter  we  idealize  soil  grains  as  two-dimensional  rods. 
This  convenient  assumption  obviously  departs  from  reality.  However,  it  is  sustainable  as 
long  as  it  provides  us  with  some  useful  hints  about  material  response.  Some  researchers 
have  already  moved  to  3D  geometry,  and  are  getting  new  insights  into  material  behavior 
(e.g.,  Thornton  and  Sun,  1994).  2D  models  are  educational  for  a  first  hands-on  experience 
of  granular  mechanics.  Many  2D  concepts  can  readily  be  extended  to  3D  modeling. 

Particles 

As  shown  in  Fig.  6,  the  set  B  represents  a  generic  assembly  of  N  rigid  particles  with 
nonlinear  interaction  at  contacts,  i.e.,  B  =  {1,...,  TV} .  The  particles,  which  are  subjected  to 

external  forces  and  moments  excluding  body  forces,  belong  to  the  set  BE  =  {N,  , 

and  the  remaining  ones  belong  to  the  set  £/.  The  sets  Bi  and  BE  are  complementary,  i.e.. 
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Figure  5.  Attributes  of  the  various  classes  of  discrete  element  methods  (adapted  from 
Cundall,  1989). 


The  particles  are  assumed  rigid  (i.e.,  they  cannot  deform).  The  assumption  of  rigid 
particles  will  be  discussed  later.  The  conservation  of  mass  implies  that  B  is  constant,  i.e., 
no  rigid  particle  should  be  lost.  Some  particles  of  Bi  may  become  part  of  Be 
forces/displacement.  However  5/  and  Be  can  vary,  provided  that  they  remain 
complementary. 
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Figure  6.  An  assembly  of  discrete  particles. 

Cylindrical  particles  are  defined  by  their  radius  and  mass  per  unit  area.  Elliptical  particles 
have  an  additional  property  to  characterize  their  aspect  ratio.  The  positions  of  rigid 
particles,  irrespective  of  their  shapes,  are  defined  by  the  center  position  and  their  rotation 
0.  The  motion  of  set  B  is  therefore  characterized  by  3 N  independent  variables. 

The  kinematics  of  the  granular  media  is  completely  defined  by  a  finite  number  of  degrees 
of  freedom.  There  is  no  need  to  use  continuous  interpolation  functions  to  represent  the 
kinematics,  unlike  in  finite  elements.  Indeed,  a  continuous  interpolation  between  discrete 
points  would  be  inappropriate  for  granular  media  because  the  displacement  is 
discontinuous  at  the  particle  contacts. 

Contacts 

Contact  mechanics  (e.g.,  Johnson,  1985)  is  a  vast  subject,  which  is  beyond  the  scope  of 
this  chapter.  We  will  only  introduce  the  basic  concepts  required  in  DEM.  With  the 
exception  of  body  forces  which  are  applied  at  the  particle  contacts,  the  external  and 
internal  forces  acting  on  particles  are  assumed  to  be  applied  at  points  and  not  to  be 
distributed  on  surfaces.  There  are  M  contact  points  belonging  to  the  contact  set  C,  i.e.,  C= 
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{1, ...,  M}.  These  contacts  can  be  subdivided  between  the  set  1, 1=  {1, Mi},  which 
corresponds  to  contact  between  particles  of  B,  and  the  set  E,  1=  {  M/+  M) ,  which 
describes  the  points  of  application  of  external  concentrated  actions.  The  sets  /  and  E  are 
complementary,  i.e.,  C  —  I U  E  =  {1,...,  M}  .  By  definition,  the  sets  Ia  and  Ea  represent 
the  contact  points  of  /and  E  on  particle  a ,  respectively.  The  subsets  Ia  and  Ea  have  the 
following  properties: 

I  =  [jla  and  E=\jEa  (1.1) 

aeB  aeB 

and 


EanEb=0  and  Ianlb={c}  for  V  a±b^B  (1.2) 

Equation  2  implies  that  two  different  sets  Ia  have  at  most  one  point  in  common.  The 
number  of  contacts  may  vary  during  a  deformation  process,  i.  e.,  /  varies  depending  on 
the  state  of  B. 


Contact  detection 

In  the  case  of  spherical  or  cylindrical  particles, 
b  when  the  following  criterion  is  met: 


there  is  a  contact  between  particles  a  and 


(1.3) 


where  x°  and  x-  are  the  center  position  of  particles  a  and  b,  and  Ra  and  Rb  their  radii, 

respectively.  The  contact  detection  criterion  becomes  more  complicated  for  elliptical 
(e.g..  Ting,  1992)  and  polygonal  (e.g.,  Cundall,  1980, 1988)  particles. 


For  large  numbers  of  particles,  the  detection  of  contacts  becomes  a  serious  computational 
issue  in  granular  mechanics.  The  most  naive  method  for  detecting  contacts  consists  of 
applying  Eq.  3  N-\  times  for  each  ball  of  an  assembly  of  AT  particles.  There  are  therefore 
Nx(N-l)  searches,  a  task  which  may  rapidly  consume  a  large  fraction  of  the  calculation 
time.  For  instance,  the  computer  would  make  10°  searches  for  10  particles,  but  10 
searches  for  106  particles,  which  is  an  excessively  large  number  of  calculations.  Several 
detection  algorithms  have  been  proposed.  The  most  efficient  ones  are  N  logN  (e.g., 
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Mujinza  et  al.,  1993).  We  will  present  one  of  the  most  commonly  used  algorithms  (e.g., 
Cundall  and  Strack,  1979). 

As  shown  in  Fig.  7,  the  space  is  divided  in  a  uniform  grid  forming  square  boxes.  The  box 
size  is  larger  than  the  diameter  of  the  largest  particles,  so  that  no  particle  covers  more 
than  four  boxes  at  once.  Assuming  that  there  is  a  maximum  number  of  M  particles  per 
box,  the  maximum  number  of  searches  is  ANM  as  a  single  particle  may  cover  4  boxes  at 
once.  This  search  technique  requires  associating  boxes  and  particles,  and  updating  this 
association  when  particles  move  across  the  grid. 

This  simple  search  technique  becomes  less  efficient  when  M  increases,  which  is  the  case 
when  the  particles  have  a  wide  range  of  sizes.  In  these  cases,  one  may  consider  alternate 
contact  algorithms. 


Contact  geometry 

The  contact  geometry  of  rigid  cylindrical  particles  is  defined  as  shown  in  Fig.  8.  In 
theory,  for  rigid  particles  which  do  not  overlap,  the  contact  is  reduced  to  the  point  of 
tangency  between  two  particles  (Fig.  8a). 


R, 


\r  = - — 

c  R,+R, 


-*A  + 


R„ 


Ra+Rb 


(1.4) 


For  particles  which  overlap  slightly,  the  contact  geometry  is  no  longer  a  point  but 
becomes  a  surface.  However,  this  contact  area  will  be  reduced  to  a  point  (Fig.  8b) 
defined  as  follows: 


1  R A 


xr  = - 

c  2 


■Rb' 


AB 


J 


,  (\  ra-rb ' 

^2  AB  j 


(1.5) 


For  rigid  cylindrical  particles,  the  contact  point  becomes  the  tangency  point  (i.e.,  Eqs.  4 
and  5  coincide).  In  the  case  of  elliptical  or  polygonal  particles,  the  contact  point  is  more 
difficult  to  define  as  described  for  elliptical  (Ting,  1992)  and  polygonal  particles 
(Cundall,  1988). 
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Figure  8.  Contact  point  between  two  cylindrical  particles  without  and  with  overlap. 
Contact  kinematics 

The  kinematics  of  the  contacts  characterizes  the  relative  motion  between  the  particles.  As 
shown  in  Fig.  9,  the  relative  displacement  of  particles  a  and  b  at  contact  point  c  is: 


A <  =  wf  - u*  +  eijk ((ob/kc  - 0)°rkac ) 


where  u°  is  the  displacement  vector  of  particle  center,  to"  the  rotation  of  particles  located 
at  x° ,  and 


,  and 


(1.7) 
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The  relative  rotation  of  particles  a  and  b: 


A co.  =(Q.  -co¬ 

in  the  case  of  cylindrical  particles, 


(1.8) 


r°c  =  *X- 


be 


=  -*X, 


CO-  =coanC3 *  and  03 i  -(0bn 3 


where  nc  is  the  unit  vector  normal  to  the  contact  area: 


(1.9) 


(1.10) 


Equation  6  becomes: 


Awf  =  wf  -  Mj°  -  «2  +  (0b^b) 

A uc2  =ub2  -u°  +ncx{coaRa  +cobRb ) 


Figure  9.  Relative  displacement  between  two  particles. 


Equation  6  describes  exactly  the  relative  displacement  between  two  particles  which  do 
not  overlap  (i.e.,  AB  =  Ra  +  Rb )•  However,  Eq.  6  becomes  approximate  in  case  of  overlap 
(i.e.,  AB<  Ra+  Rb).  Fortunately,  this  overlap  is  usually  very  small  in  most  cases,  and  is 
thought  to  have  no  significant  consequences. 


Contact  actions 

The  actions  at  contact  c  are  represented  by  contact  forces  f°c  and  contact  moments  m°c . 
Both  forces  and  moments  are  applied  at  the  contact  points.  In  many  instances,  the  contact 
moments  are  neglected,  and  the  contact  actions  are  reduced  to  forces  alone.  This 
assumption  may  apply  to  small  contact  areas  between  cylindrical  and  spherical  particles, 


which  cannot  transmit  significant  moments  due  to  their  small  size.  However,  this 
assumption  may  become  questionable  for  particles  of  arbitrary  shapes  and  large  normal 
contact  forces. 


Contact  relations 

Various  types  of  contact  relations  between  particles  were  recently  reviewed  in  (Misra, 
1995),  including  contacts  between  smooth,  spherical,  non-spherical,  cylindrical,  and  non- 
cylindrical  elastic  particles  with  friction  and  surface  adhesion,  rough  elastic  particles,  and 
viscous  bridge.  Hereafter,  we  will  only  review  the  normal  stiffness  between  spherical  and 
cylindrical  particles,  and  some  findings  on  contacts  with  friction. 


Elastic  contact  between  smooth  particles 

The  distribution  of  contact  pressure  proposed  by  Hertz  (1882)]  is: 

P(r)  =  p0<Jl-{r/df  O-12) 

where  po  is  the  maximum  contact  pressure,  a  is  the  radius  of  the  circular  contact  area,  and 
r  is  the  polar  coordinate.  The  total  load  P  is  related  to  the  contact  pressure  through: 

P  =  p(r)2ttrdr  =  ^p0Jta2  (1-13) 

Therefore  the  maximum  pressure  po  is  3/2  times  the  mean  pressure  pm.  The  radius  of  the 
contact  area  is: 

a  =  ^0JR/2£*=3j|p  (1-14) 

where: 
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(1.15) 


In  the  case  of  identical  elastic  properties,  E*=E/2(l-v*).  The  mutual  approach  of  distant 
points  in  the  two  solids  is: 


«-«p0/22r=^--J-^=-  «  r  =  ±E-R2  (8/Rf 


R  VI 6RE 


The  secant  normal  stiffness  kn  is  therefore  load-dependent: 


(1.16) 


*.  =4  =!&/>« 
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In  the  case  of  cylindrical  solids,  the  half-width  of  the  contact  area  is: 


(1.17) 


(1.18) 


where  P  is  now  a  force  per  unit  length  (i.e.  [P]  =  F  L'1).  The  maximum  contact  pressure 
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The  relative  displacement  8 is  given  in  (Johnson,  1985): 


(1.19) 


8  =  -^-(LnU7tRE'/p)-l) 
nE 

The  secant  normal  stiffness  k„  is  therefore  load-dependent: 


(1.20) 


P  nE 

”  “  8  ~  LnUnRE'/P)-l 


(1.21) 


Figure  10  shows  the  variation  of  normal  load  P  with  normal  relative  displacement  8  in 
the  case  of  spherical  and  cylindrical  contact.  Figure  1 1  shows  the  corresponding  secant 
normal  stiffness  k„.  The  relative  displacement  8  is  normalized  with  the  particle  radius  R. 
The  normal  load  P  is  normalized  with  RE*  for  cylindrical  particles  and  RE*2  for 
spherical  particles.  In  the  two-dimensional  case  (i.e.,  cylinder),  [P]=  FL'1,  and  in  the 
three-dimensional  case  (i.e.,  sphere)  [P]=  F. 
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Figure  10.  Variation  of  normal  load  with  relative  displacement  <VR  for  cylindrical  (2D) 
and  spherical  (3D)  particles 
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Figure  11.  Variation  of  normal  secant  stiffness  corresponding  to  Fig.  10. 

Table  1  lists  typical  values  of  elastic  properties  for  various  particle  materials.  These 
values  are  useful  to  define  realistic  values  of  contact  stiffness  between  particles. 


001  0.01  0.1  1 
SIR 
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Table  1.  Elastic  properties  for  particle  materials  (after  (Bardet,  1997;  Johnson,  1985). 


Material 

Young's  modulus  (Gpa) 

Poisson  ratio 

Perpex 

3 

0.38 

Glass 

55 

0.25 

Steel 

200 

0.28  -  0.29 

Aluminium 

55-76 

0.34  -  0.36 

Duraluminium 

74 

0.32 

Cast  Iron 

113 

0.25 

Tungsten  Carbide 

732 

0.22 

Amphibolite 

93  - 121 

0.28  -  0.30 

Anhydrite 

68 

0.30 

Diabase 

87-117 

0.27  -  0.30 

Diorite 

75  - 108 

0.26  -  0.29 

Dolomite 

110-121 

0.30 

Dunite 

149-183 

0.26  -  0.28 

Feldspathic  Gneiss 

83-118 

0.15-0.20 

Gabbro 

89-127 

0.27-0.31 

Granite 

73-86 

0.23  -  0.27 

Limestone 

87-108 

0.27  -  0.30 

Marble 

87-108 

0.27  -  0.30 

Mica  Schist 

79-101 

0.15-0.20 

Obsidian 

65-80 

0.12-0.18 

Oligoclasite 

80-85 

0.29 

Quartzite 

82-97 

0.12-0.15 

Rock  salt 

35 

0.25 

Slate 

79-112 

0.15-0.20 

Ice 

7.1 

0.36 

Deformable  and  rigid  particles 

As  previously  mentioned,  soil  particles  were  assumed  to  be  rigid  and  the  contact 
deformable.  However,  in  reality,  soil  particles  are  not  rigid;  they  deform  when  they  are 


subjected  to  contact  forces.  The  assumption  of  rigid  particles  and  deformable  contacts  is 
acceptable  as  long  as  the  contact  deformations  represent  the  particle  deformations.  This 
condition  may  be  met  for  particles  with  simple  geometry  (e.g.,  spheres  and  cylinders) 
undergoing  elastic  deformations.  However,  this  condition  may  not  be  met  for  complex 
deformation  patterns  and  inelastic  deformation,  and  distributed  contact  loads  (e.  g., 
Cundall,  1989).  The  case  of  deformable  particles  is  not  considered  hereafter.  One  may 
refer  to  Cundall  (1980)  and  Shi  (1982)  to  account  for  elastic  deformation  of  particles. 

We  will  illustrate  the  fact  that  the  contact  deformation  may  represent  the  particle 
deformation  in  some  simple  cases.  This  will  be  demonstrated  by  considering  the 
compression  of  the  cylinder  of  Fig.  12,  which  is  subjected  to  diametrically  opposed 
concentrated  forces. 

The  stress  distribution  in  the  cylinder  is  given  in  Timoshenko  and  Goodier  (1951).  The 
stresses  at  point  A  are: 


P  1  2{a]  +2z2)  Az  P  12  2 

- L  +  —  and  a,  = - - ,  ,  ■— 

K  R  a* Jaj2  +z2  ai  J  JtyR  2R-z  yjoi+z 


(1.22) 


In  plane  strain,  the  vertical  strain  is: 


(1.23) 
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Figure  12.  Compression  of  a  cylinder  due  to  diametrically  opposed  concentrated  loads. 

The  compression  of  the  upper  half  of  the  cylinder  0,C  is  found  by  integrating  &  from  z  = 
0  to  z  =  R,  where  a  «R,  to  give: 


«,  =  f’e.dt  =  f(1  V|’)(2to(4R/<.l  )-l) 

J0  7tt 


(1.24) 


A  similar  expression  is  obtained  for  the  compression  of  the  lower  half  of  the  cylinder  so 
that  the  total  compression  of  the  diameter  through  the  mid-points  of  the  contact  areas 
O1O2  is: 


8  =  2P(l  V  l(Ln(4R/ai )+  Ln{4R/a2 )- 1) 
7lE 


(1.25) 


Equation  25  can  be  used  to  derive  Eq.  20.  As  clearly  indicated  by  Eqs.  22  to  25,  the 
contact  stiffness  originates  from  the  particle  deformation.  In  this  particular  case,  the 
contact  deformation  accounts  for  the  particle  deformation. 

Frictional  contact 

In  soil  mechanics,  the  frictional  characteristics  between  two  particles  is  commonly 
thought  to  be  a  basic  material  property,  which  can  easily  be  determined  from  laboratory 
experiments  and  tables  of  physical  constants.  However,  this  common  belief  is 
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unfortunately  unfounded.  The  determination  of  friction  between  two  particles  is  still  a 
complicated  problem  as  described  in  Singer  and  Pollock  (1992).  As  shown  in  Table  2,  the 
values  of  friction  angle  vary  not  only  with  mineral  type  but  also  on  the  contact 

cleanliness,  water  content,  and  level  of  normal  load. 

Contact  models 

Several  relations  were  proposed  for  relating  the  contact  actions  and  contact  kinematics. 
We  will  review  only  two  basic  ones:  the  elastic-perfectly  plastic  relations  with  and 
without  rotational  stiffness  and  friction. 

Elastic-perfectly  plastic  contact 

Figure  13  represents  the  force-displacement  relationship  between  two  cylindrical 
particles  which  is  used  to  simulate  the  intergranular  behavior.  The  contact  relationship  is 
activated  when  two  disks  overlap.  The  contact  geometry  between  two  disks  is 
characterized  by  a  contact  point,  and  a  contact  direction  which  passes  through  the  centers 
of  the  particles  in  contact.  N  and  S  denote  the  projections  of  the  contact  force  that  are 
respectively  tangential  and  normal  to  the  contact  direction.  The  value  of  the  contact  force 
at  time  t+ At  is  calculated  from  its  value  at  time  t  by  using  the  following  relation 


N(t  +  At)  =  N(t)  +  k„An 
S(t  +  At)  =  S(t)  +  ksAs 


(1.26) 


where  k„  and  ks  are  the  tangential  normal  and  tangential  stiffness,  respectively,  and  An 
and  As  are  the  normal  and  tangential  components  of  the  relative  displacement  of  the 
contact  between  time  t  and  t+At.  In  the  case  of  linear  elastic  contacts,  k„  and  ks  are 
constant.  In  the  case  of  Hertz  contact,  k„  and  ks  vary  with  the  contact  load  N  and  S.  In 
general  kn  and  ks  are  assumed  to  be  equal.  More  realistic  expressions  are  reported  in 
Misra  (1995). 
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Table  2.  Interparticle  friction  angle  fa  measured  from  laboratory  tests  (after  Misra,  1995). 


Material 

Type  of  test 

Conditions 

Value  <pu 

Reference 

Biotite 

Along  cleavage  faces 

Dry 

14.6-17.2 

Horn  and  Deere, 

1962 

Saturated 

7.4 

— 

Calcite 

Block  on  block 

Dry 

8 

_ 

Water-saturated 

34.2 

— 

Chlorite 

Along  cleavage  faces 

Dry 

19.3-27.9 

— - 

Saturated 

12.4 

— 

Feldspar 

Block  on  block 

Dry 

6.8 

— 

Water-saturated 

37.6 

— 

Direct  shear  box,  free  particles  on  flat  surface 

36 

El-Sohby,  1969 

Free  particle 

37 

Lee,  1966 

Particle-plane 

Saturated 

28.9 

Procter  and 

Barton,  1974 

Feldspar 

Direct  shear  box,  fixed 

Dry 

6 

Home,  1965 

(microline) 

particles  on  flat  surface 

Feldspar  (microline) 

Water-saturated 

37 

— 

Glass 

Free  particles  on  plate 

Dry 

10-12 

Gray,  1960 

Water-saturated 

17 

— 

Three  balls  on  glass  plate 

Dry,  Low  load 

9 

— 

Water-saturated  and  low  load 

19 

— 

After  Redrying 

16 

— 

Glass  ballotini 

Ball  on  ball 

Dry,  low  load 

2 

Skinner,  1969 

Dry,  high  load 

4 

— 

Flooded,  low  load 

28-40 

— 

Flooded,  high  load 

38-40 

— 

Dry,  low  load 

3 

— 

Dry,  high  load 

7 

— 

Direct  shear  box,  free 

Water-saturated 

17 

Rowe,  1962 

particles  on  flat  surface 

Dry,  tested  in  dry  nitrogen 

9 

Tong,  1970 

Acetone  cleaned,  tested  in  dry 

16 

— 

nitrogen 

Trichloroethylene,  acetone,  detergent 

21 

— 

rinses 

Water-saturated 

15 

— 

Cleaned  with  soap,  water  and  acetone 

15 

— 

Particles  fixed  with  wax  after  initial  sliding 

14-15 

— 

Muscovite 

Along  cleavage  faces 

Dry 

16.7-23.3 

Horn  and  Deere, 

1962 

Saturated 

13 

— 

Phlogopite 

Along  cleavage  faces 

Dry 

14-17.2 

— 

Saturated 

8.5 

— 

Phosphor-bronze 

Ball  on  ball 

Water-saturated 

21 

Parikh,  1967 
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Material 

Quartz 


Type  of  test 
Block  on  block 


Reference 
Bromwell,  1966 
Horn  and  Deere,  1962 


Block  over  particle  set  in 
mortar 

Direct  shear  box  fixed 
particles  on  flat  surface 


Particle-particle 

Particle-plane 

Particles  on  polished  block 
Three  fixed  particles  over 
block 

Quartz  (clean)  Direct  shear  box,  fixed 
particles  on  flat  surface 


Quartz  (milky) 


_ Conditions _ 

Dependent  on  surface  condition 
Dry 

Water-saturated 

Dry 

Moist 

Water-saturated 

Dry 

Moist  and  water-  saturated 
Dry 

Water-saturated 
High  load 
Low  load 
Air  dried 
Atmospheric 
Water-saturated 
High  vacuum 
Water-saturated 
Moist  and  water-saturated 
Dry,  tested  in  dry  nitrogen 
Saturated 
Saturated 
Dry 

Water-saturated 

Water-saturated 

Dry 


Value  i 
0-45 
7.4 

24.2 

6  Tschebotarioff  and 

Welch,  1948 
24.25  — 

24.25  — 

6  — 

25  — 

11  Penman,  1953 

33  — 

19  — 

29  — 

22  Bishop,  1954 

28  Brogliato,  1996 

23-30  [175] 

38  Bromwell,  1966 

26  Tong,  1970 

28  — 

15  — 

26  Procter  and  Barton,  1974 

22.2  — 

7.4  — 

22-31  Rowe,  1962 

21-27  Hafiz,  1950 

6  Horn  and  Deere,  1962 


Water-saturated 

23 

— 

Amylaraine 

31 

— 

Carbontetrachloride 

11 

— 

Dry 

9 

— * 

Water-saturated 

27 

— 

Quartz  (rose) 

Dry 

Water-saturated 

7 

24 

_ 

Steel 

Ball  on  ball 

Dry,  polished  and  cleaned  with 
carbontetrachloride 

7 

Gray,  1960 

Direct  shear  box,  free 
particles  on  flat  surface 

Air 

Water-saturated 

7 

9 

Rowe,  1962 

Free  particles  on  plate 

Water-saturated,  dry  and  cleaned  with 
carbontetrachloride 

8.8 

Gray,  1960 

Light  load  apparatus 

Dry,  polished 

Dry,  polished  and  cleaned  with 
carbontetrachloride 

9 

9.5-12 

— 

Rods  on  rods 

Dry,  cleaned  with  carbontetrachloride 

9-14 

— 

Steel  balls 

Friction  apparatus 

Dry 

16-32 

Skinner,  1969 

Zircon 

Direct  shear  box,  free 
narticles  on  flat  surface 

Water-saturated 

23 

El-Sohby,  1969 
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The  contact  force  of  Fig.  13  obeys  the  Coulomb  friction  law: 


S(t)  <  tan0MW(f)  +  c  (1-27) 

where  and  c  are  the  friction  angle  and  cohesion  between  two  disks,  respectively. 


Figure  13.  Elastic  perfectly-plastic  model  for  contact  without  rotational  stiffness  and 
friction. 

Linear-elastic  perfectly  plastic  contact  with  rolling  stiffness  and  friction 

In  idealized  granular  materials,  rolling  friction  is  included  at  the  particle  contacts  by 
generalizing  the  elastic  perfectly-plastic  contact  relation.  Such  a  generalized  model  is 
schematized  in  Fig.  14,  including  its  normal,  tangential,  and  rotational  stiffnesses,  and 
rolling  and  sliding  frictions.  At  the  particle  contacts,  the  increments  of  normal  force, 
shear  force,  and  moment  are: 

AFn  =  knAn,  AFs=ksAs,  and  AM  =  RkeA6  (1-28) 

where  kn,  ks,  and  ke  are  the  normal,  tangential,  and  rotational  stiffness  of  grain  contacts, 
respectively.  R  is  the  average  particle  radius,  An  and  As  are  the  normal  and  tangential 
relative  contact  displacement,  and  AO  is  the  relative  rotation  of  disks,  which  is  the 
rotation  of  the  contact  surface.  Eq.  28  holds  provided  that  the  normal  and  tangential 
forces,  and  the  moment  at  contact,  obey  the  following  generalized  Coulomb  friction  law: 


b  =  -R2 \A62  -  fi\ sign  {(A0,  -  j3)(A02  -  p)} 


(1.29) 


# 


+ 


♦ 


where  is  the  sliding  friction  angle,  and  the  rolling  friction  angle.  The  rotational 
stiffness  kg  is  based  on  an  analytical  expression  derived  from  the  two-dimensional  theory 
of  elastic  contact  (Meftah  et  al.,  1993).  ke  increases  proportionally  to  the  cylinder  radius 
and  the  normal  load,  which  is  in  agreement  with  experimental  results  on  hard  rubber 
cylinders  (Bardet  and  Huang,  1993). 

Governing  equations  of  statics 

The  equilibrium  of  forces  on  particle  a  is: 

!/“  +  £/“=  0  *  =  1.2,3  (1.30) 

c€/„  eeEa 

where  /« is  the  internal  force  at  contact  c,  and  fae  the  external  force  at  point  e.  The 
equilibrium  of  moments  about  center  of  particle  a  is: 

I(<+^'-r/r)+i(<+^'T/i”)=o  <=1.2,3  (1.3D 

where  m.c is  the  internal  moment  at  contact  point  c;  mf  the  external  moment  at  point  e\ 
eijk  the  alternating  tensor;  rf  =xcj-xaj,  rf  =  x]  -  xbj  ;  x°  and  x)  the  center  coordinates 
of  particles  a  and  b;  and  xc}  the  position  of  contact  point  c. 


# 
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# 


# 


# 


» 


Figure  14.  Representation  of  generalized  contact  with  rotational  stiffness  and  rolling 
friction. 

Boundary  conditions 

There  are  several  types  of  boundary  conditions,  which  can  be  specified  on  an  assembly  of 
particles,  including  (1)  prescribed  displacement,  velocity,  acceleration,  and  force 
boundary,  (2)  periodic  boundary,  (3)  rigid  boundary,  and  (4)  flexible  boundary. 

Prescribed  force/displacement 

As  in  the  boundary  value  problems  of  continuum  mechanics,  the  boundary  conditions  in 
granular  mechanics  can  be  either  prescribed  displacement/  velocity/  acceleration,  or 
prescribed  force.  External  forces  and  moments  can  be  applied  to  any  point  of  particles. 

The  motion  of  a  group  of  particles  can  conveniently  be  prescribed  by  a  cluster.  Clustered 
particles  are  subsets  of  particles  that  move  as  a  rigid  body  and  are  useful  to  represent 
rigid  objects.  As  shown  in  Fig.  1 5,  the  particles  may  overlap,  and  move  as  a  solid  object. 
The  motion  of  the  ath  particle  in  a  cluster  is  defined  by  its  translation  and  rotation  AQa 

as  follows: 

u°  =  U,  +  eljkA6  (X,  -  x° )  and  AQa  =  A9  (1 .32) 

where  U,  is  the  translation  vector  and  A 6  is  the  rotation  of  a  reference  point  for  the 
cluster. 


6 
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Figure  1 5.  Rigid  cluster  of  particles. 

Periodic  boundaries 

Periodic  boundaries  are  illustrated  in  Fig.  16.  The  particles  leaving  segment  AB  are 
reintroduced  on  segment  CD.  The  periodic  boundary  conditions  can  be  used  in  both 
horizontal  and  vertical  directions,  therefore  filling  the  complete  space  with  particles. 
Periodic  boundary  conditions  have  extensively  been  used  for  determining  average 
continuum  quantities  (e.g.,  stress  and  strain)  free  of  the  heterogeneities  caused  by  rigid 
boundaries.  However,  periodic  boundaries  introduce  kinematic  constraints  in  some 
circumstances.  These  conditions  are  met  when  the  deformation  patterns  have  a  length 
scale  that  is  a  sub-period  of  the  box.  One  of  these  adverse  effects  is  observed  for  shear 
bands  as  explained  in  Bardet  and  Proubet  (1991). 


Spatial  period 


B  D 


Figure  16.  Periodic  boundaries. 

Rigid  walls 

Rigid  walls  simulate  the  loading  platens  which  transmit  loads  to  laboratory  samples. 
They  are  usually  made  of  a  single  segment,  but  can  also  be  made  of  several  connected 


* 


* 


segments  for  generating  various  polygonal  shapes,  such  as  superficial  footings. 
Displacement,  force,  or  moment  can  be  specified  for  rigid  walls. 

Flexible  membranes 

Flexible  membranes  were  introduced  in  Bardet  and  Proubet  (1991)  to  simulate  the 
flexible  membranes  used  in  the  triaxial  test.  Forces  and  moments  are  externally  applied  to 
the  boundary  particles  by  specifying  the  prescribed  stress  tensor.  The  forces  distributed 
on  the  boundary  are  calculated  from  the  unit  vectors  normal  to  the  boundary  segments 
and  the  prescribed  stress  tensors.  For  instance  in  Fig.  17,  the  force  F°  applied  to  the 

particle  center  O  is: 

F°  =  BCovnBC  ( OC  +  CBI2)I0A  +  CDa^n™  +  DEoynEE  ( OD  +  DE/2)/OF (2.33) 

where  a y  is  the  prescribed  stress  tensor  and  nBC ,  nCD  ,  and  „DE  are  the  unit  vectors  normal 

to  segments  BC,  CD  and  DE.  The  flexible  membrane  can  sustain  large  deformation. 
Particles  are  inserted  into  the  flexible  boundary  chain  as  they  attempt  to  force  their  way 
between  two  particles  of  the  flexible  boundary. 


o 


Figure  17.  Flexible  membrane  for  stress-controlled  boundaries. 
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PART  II.  TRANSITION  FROM  DISCRETE  TO  CONTINUOUS  MEDIA 


Background 

The  definition  of  stress  is  one  of  the  critical  problems  for  understanding  material 
instability  and  the  assumptions  of  higher-order  continuum  theories.  The  definition  of 
stress  in  granular  media  is  a  controversial  topic  in  Mechanics.  Some  researchers  (e.g., 
Bogdanova-Bontcheva  and  Lippmann,  1975;  Chang  and  Ma,  1991;  Kanatani,  1979;  and 
Mulhaus  and  Vardoulakis,  1987)  claim  that  stress  tensor  is  not  symmetric  in  granular 
media,  and  that  couple  stresses  are  important  to  understand  material  instability  such  as 
shear  banding.  Others  (e.g.,  Christoffersen  et  al,  1981;  and  Cundall  and  Strack,  1978) 
affirm  that  the  stress  asymmetry  is  absent  or  negligible  for  all  practical  purposes,  and 
unnecessarily  complicates  the  description  of  the  mechanical  behavior  of  granular  media. 

The  controversy  about  the  asymmetry  of  stress  and  existence  of  couple  stress  is  not 
specific  to  granular  media.  Couple  stresses  were  proposed  in  metals  and  fracture 
mechanics  to  regularize  the  stress  intensity  at  crack  tips  (e.g.,  Sternberg,  1968),  but  there 
is  not  yet  a  convincing  experimental  evidence  for  couple  stress  (e.g.,  Diepolder  et  al., 
1991). 

Part  II  re-examines  the  definition  of  stress  in  granular  materials,  and  establishes  the 
conditions  under  which  there  may  be  couple  stresses  and  asymmetric  stress.  Following 
the  introduction,  the  second  and  third  sections  review  the  basic  equations  of  granular  and 
equivalent  continuous  media.  The  third  section  re-examines  the  definition  of  stress  from 
virtual  work  and  statics.  Finally,  the  last  section  gives  an  example  illustrating  the  stress 
asymmetry  in  granular  media. 

Granular  medium 

Definition 

As  shown  in  Fig.  1,  the  volume  V  is  filled  with  N  particles,  some  of  which  are  subjected 
to  external  forces  or  moments  applied  from  the  exterior  of  volume  V.  The  particles  are 
grouped  in  the  set  5=  {1,  The  forces  and  moment  acting  on  the  particles  of  B  are 


concentrated  at  M  points  of  set  C—  (1, ... ,M}.  As  shown  in  Fig.  1 ,  the  subset  I  represents 
the  contact  points  between  two  particles  of  B,  whereas  the  set  E  denotes  the  points  where 
external  actions  are  applied: 

7  =  {1,...,M,},  E  =  {MI+l,...,M),  and  C  =  I  v  E  =  (2-1) 

The  sets  Ia,  Ea,  and  Ca  denote  the  contact  points  on  particle  a  corresponding  to  internal 
actions,  external  actions,  and  all  actions,  respectively.  The  sets  Ca,  Ia>  Ea,  I,  E,  and  C  are 
related  as  follows: 

c = (Jc„,  C,  =  I.uE„  /=|_K>  ■>"“  <2-2> 

aeB  oeB  aeB 

The  intersections  of  Ia  and  Ea  for  two  different  particles  are  either  empty  or  reduced  to  a 
single  point  c : 

EanEb=0  and  Ia  n  Ib  =  {c}  for  V  a*  be  B  (2.3) 

The  particle  assembly  is  in  equilibrium  when  each  particle  is  in  equilibrium.  The 
equilibrium  of  internal  and  external  forces  acting  on  particle  a  is: 

^f;c  =0  /  =  1,2,3  (2-4) 

ceCa 

where  f°c  is  the  force  at  contact  c.  The  equilibrium  of  moments  about  the  center  of 
particle  a  is: 

-*;>/,”)=  0  /-  1,2,3  (2.5) 

ceCa 

where  m°c  represents  the  internal  or  external  moments  at  contact  point  c;  eiJk  is  the 
permutation  symbol  used  for  vector  cross  product;  and  x°  and  xCj  are  the  coordinates  of 
particle  center  a  and  contact  point  c,  respectively. 
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Virtual  work  in  granular  media 


As  shown  in  Fig.  1,  the  kinematics  of  granular  media  is  represented  by  the  displacement 
8u°  of  the  particle  centers,  and  the  particle  rotation  86° .  After  multiplying  Eqs.  4  and  5 

by  any  virtual  displacement  8u°  and  rotation  86° ,  and  summing  for  all  the  particles  of 
volume  V,  one  obtains  the  following  relation: 

X X  (r&r + (< + %(*;  -  *;>/*")»"  )= 0  <2-6> 

aeBceC0 

After  transforming  the  double  sum  for  Ca  and  B  of  Eq.  6  into  two  separate  sums  for  /  and 
E,  and  noting  that  contact  forces  and  moments  are  opposite  at  internal  contact  (i.e., 
fc  _  j-ac  _  _  j-bc  an(j  mc  _  mac  _  _mbc  ^  one  obtains  the  principle  of  virtual  work: 

8W,d  +  8Wed=  0  (2.7) 

where  the  work  8WE  done  by  external  forces  and  moments  and  the  work  8W,D  done  by 
internal  forces  and  moments  are: 

+  and  8WED=-^{f;8u°  +  mJS0;)  (2.8) 

cel  eeE 

As  shown  in  Fig.  2,  A  8u°  and  A  86°  are  the  relative  displacement  and  rotation  of  the  two 
particles  a  and  b  at  their  contact  point  c,  respectively: 

A 8u°  =  8u>  -  8u°  +  eiJk (< 56 ) (xtc  -  xbk )-  86° (x°  -  x° ))  and  A 86°  =  86*  -  86°  (2.9) 

8u-  and  86-  are  the  displacement  and  rotation  of  the  points  e  of  application  of  external 
forces  and  moments.  The  variational  displacements  and  rotations  8u°  and  86°  can  be 
selected  arbitrarily.  In  particular,  they  can  be  chosen  as  follows: 

8u°=ai+biJx°+ciJkx°x°k  and  86°=at  +  l %x°  i=  1,2,3  (2.10) 

where  a,,  by,  Cyk,  Ok,  and  fiy  are  arbitrary  coefficients.  By  using  Eq.  10,  A8u° ,  A 86°  and 
8u‘  become: 
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A 5m, c  =  bij{xbj-x°)+cijk(xbjxbk  -XjXak)-cCjeiJk(xbk  -xak)+  Pj,eijk(xb(xck -xbk)-xf{xck  -x,°))(2.1 1) 


A80‘  =  Pij{xb-x°) 


(2.12) 


8u‘  =  Su°  +  etjk8e°  {xk  -  x? )  =  a,  +  bijXa;  +  c,jkx°exake  +  eijka]  -  xake )+  el]k pjtx?  (xek  -  xake  )(2. 1 3) 

where  xae  corresponds  to  the  center  of  particle  a  where  contact  e  takes  place.  By  using 
Eqs.  11  to  13,  8WjD  and  8Wk  become: 


8W,D  = 


8W. 


D  „ 


bj'Lf.ixj -x°i)+cijk'Lfi(xbA -XjX^-a^eyJ^xh -xak) 

C€/  /  /  ee/  \  ee/t  \  (2.14) 

+  P jl^ipijkff  (X/  (X*  ~xk)~xl  (Xk  -  Xk  ))+  mj(Xl  ~XlV 

cel 

f'xr -oc^e, „/;<.< -*n+™;) 

eeE  eeE  eeE  eeE  (2.15) 


eeE 


Because  Eqs.  7, 14  and  15  hold  for  arbitrary  values  of  ah  by,  cyk  at  and  py,  the  following 
relations  are  obtained: 


£/‘=0  i=  1,2,3 

eeE 

(2.16) 

cel  eeE 

i,j=  1,2,3 

(2.17) 

=-XM" 

cel  eeE 

/=  1,2,3 

(2.18) 

cel 

-xZ))+m;(x‘-x")h'ZMrx“  i,j~  1,2,3 

eeE 

(2.19) 

'LfttfA  -*w=1fr= 

cel  eeE 

ca/xake  i,  j,  k=  1,2,3 

(2.20) 

where  M°e  is  the  external  moment  acting  on  particle  a  about  the  center  of  particle  a: 


M?  =e,]k{x)-xa;)rk+m‘  (2.21) 

Equation  16  translates  the  equilibrium  of  external  forces  applied  to  the  whole  assembly  of 
particles.  By  using  Eqs.  17  and  18,  one  can  derive  the  equilibrium  of  external  moments 
about  the  coordinate  origin  for  the  assembly  of  particles,  i.e.: 
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(2.22) 


£femX/„e +<)= o 

ee£ 

Therefore,  Eq.  18  becomes: 

-*;>/;  -  -s«r  - 1  «-,</; 

ce/  ee£  ee£ 

For  a  volume  Fto  be  in  equilibrium,  the  sum  of  external  moments  about  a  common  point 
must  vanish  (i.e.,  Eq.  22).  However,  it  is  emphasized  that  the  sum  of  external  moments 
about  different  particle  centers  (i.e.,  )  is  not  necessarily  equal  to  zero.  M°e  results 

eeE 

from  not  only  contact  moments  m‘  but  also  contact  forces  f-  .  It  may  be  different  from 
zero  even  where  there  is  no  contact  moment. 

Continuum  for  granular  media 

In  the  continuum  equivalent  for  granular  media,  the  traction  vector  7}  and  moment  vector 
m,  acting  on  the  unit  surface  of  unit  normal  vector  n,  is  related  to  the  Cauchy  stress  tensor 
Gij  and  the  couple  stress  tensor  through: 

Ti  =  Oj,nj  and  w,  =/xy/«7  i=  1,2,3  (2.24) 

In  the  absence  of  external  body  force  and  moment  per  unit  volume,  the  equations  for 
equilibrium  of  internal  stress  and  couple  stress  are: 

(Jjij  =  0  i  =1,2,3  (2.25) 

Pjij  +  e,ki°ki  =  0  i=  1  (2.26) 

The  kinematics  of  the  equivalent  continuum  is  defined  by  the  fields  of  displacement 
vector  Suj  and  rotation  86i,  which  describe  the  motion  of  particle  centers,  i.e.: 

8ul(x°)  =  8u°,  and  86^°)=  86°  VaeB  (2.27) 

By  multiplying  Eqs.  25  and  26  by  any  variational  fields  8uj  and  86,,  and  integrating  over 
the  volume  V,  one  obtains  the  following  relation: 


(2.28) 


+elklald)8ei)=  0 

By  invoking  Gauss  theorem,  the  principle  of  virtual  work  is  obtained: 

SW,+8We=  0  (2.29) 

where  the  virtual  work  SWe  of  external  forces  and  moments  and  the  virtual  work  SJF/of 
internal  stresses  are: 

dWj  =  l (a fi (SU'  j  +  euk86k )+  n/O.j )iV  and  8WE=-\s{Tl8u,  +  ml86)iS  (2.30) 
By  choosing  and  8Qi  as  specified  in  Eq.  1 0,  Eq.  30  becomes: 

8We  =  -a,J5  TtdS  -  btJ\s  T'XjdS  -  ciJk\s  T,XjXkdS  -  a,  £  mtdS  -  mlXjdS  (2.31) 

8W}  =  b^Cj.dV  +  c^(<V*  +okixJ)dV-ocievkl<TjkdV  +  Pyj^+e^jO^V  (2.32) 
Because  Eqs.  29,  31  and  32  hold  for  any  values  of  a„  by,  CyhOCj,  and  A/,  the  following 


relations  are  obtained: 

[7^5  =  0  /=  1,2,3  (2-33) 

J  s' 

J/ TijdV  =  \sxiTJdS  i,j=  1,2,3  (2.34) 

e0k  l °jkdV  =  -£ mtdS  i  =  1 ,2,3  (2.35) 

l (A,  +  eJk,Xi°Ik  )dV  =  £ x^dS  ij  =  1,2,3  (2.36) 

+  OkiXj)dV  =  jT^dS  i, j  =  1 ,2,3  (2.37) 


Equation  33  implies  the  equilibrium  of  external  forces.  Equation  34  represents  the 
average  stress  otj  in  volume  V: 

a  =—\  audV=-\  XiT.dS 
‘J  y  Jv  ,J  yJs  '  J 

Equation  35  is  useful  to  examine  the  symmetry  oi  a y\ 


(2.38) 


eikj°kJ=-^\sm,ds  /=1’2’3  (2-39) 

elkjakj  =  0  when  the  stress  tensor  is  symmetric  (i.e.,  atJ  =  ajt ).  However,  £  mtdS  is  not 

necessarily  equal  to  zero  when  there  are  moments  at  the  external  boundary.  As  previously 
mentioned,  these  external  moments  in  granular  media  may  result  from  contact  forces 
without  contact  moments.  Finally,  Eq.  36  becomes: 


y  jy  (j^ij  "*■  ejklXPlk  A ij  ejkl^‘ilk 

where  %Jk  is  the  average  moment  of  stress: 

^‘jk  =  Jpfy  XPjkdV 

In  summary,  the  internal  work  becomes: 


(2.40) 


(2.41) 


8W,  =  V(b0.ajt  +  cljk(ZkJ,  +  'Ljki)-ajetjkajk  +  A,(/T,,  +eiklZj!k))  (2.42) 

Definition  of  average  stresses  in  granular  media 

The  average  stresses  in  the  equivalent  continuum  are  defined  by  postulating  that  the 
granular  and  continuous  media  produce  identical  internal  and  external  works: 


SW,D  =  8Wj  and  8WCD  =  8WE 
Average  stress 


(2.43) 


Because  Eq.  43  applies  to  arbitrary  values  of  by,  the  average  stress  is: 


X*r/; 

*  cel  '  eeE 


(2.44) 


Equation  44  is  identical  to  those  derived  by  Weber  (1966),  Goddard  (1977), 
Christoffersen  et  al.  (1981),  and  Rosenberg  and  Selvadurai  (1981).  In  the  case  of 
spherical  and  cylindrical  particles  a  and  b,  which  are  in  contact  at  point  c  with  unit 
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normal  vector  n-  (i.e.,  xf  -  x°  =  (Ra  +  Rb)n ■  ),  due  to  the  opposite  sign  of  contact 

forces  and  contact  normals  (i.e.,  ff  =  -fcb  and  n°°  =  -nf ),  Eq.  44  becomes: 


«.k/; iw 

K  ce/ 


aeB  cel  a 


Symmetry  of  average  stress 


(2.45) 


Because  Eq.  43  applies  to  arbitrary  values  of  (Xj,  one  obtains: 


V*#  -  ~z«r  -  i=l-2-3 

^  cel  *  eeE  *  eeE 


(2.46) 


Equation  46,  which  can  also  be  obtained  directly  from  Eq.  44,  is  useful  to  determine  the 
amplitude  of  stress  asymmetry.  This  amplitude  can  also  be  characterized  by  cTy  -  tr),  as 

follows: 


^  eeE 


ae 

k 


(2.47) 


eeE 


Equation  47  implies  that  the  average  stress  may  be  asymmetric,  even  when  there  is  no 
moment  at  contacts  (i.e.,  m‘=  0).  The  asymmetry  results  from  the  sum  of  the  external 

moments  that  are  created  by  external  forces  fe  about  the  particle  centers. 


The  amplitude  of  stress  asymmetry  increases  with  the  area  S  on  which  the  external 
moments  are  applied,  but  decreases  with  the  volume  size  V.  If  the  external  moments  are 
assumed  to  have  bounded  values,  the  amplitude  of  stress  asymmetry  decreases  with  V/S. 
When  V  — »  the  effects  of  external  moments  vanish,  and  the  average  stress  is 

symmetric.  This  applies  with  or  without  contact  moments. 


In  the  case  of  spherical  and  cylindrical  particles,  Eq.  47  becomes: 


z,r”,,=  'yTLR. k/r - </rh v !«, + Ki<fr - <fr)  (248) 


ae  B  ce  la 


cel 


Equation  48  shows  that  a(>  is  not  necessarily  symmetric  when  the  particles  are  spheres  or 

cylinders  of  identical  radius.  This  result,  which  is  in  disagreement  with  Caillerie  (1991) 
and  Chang  and  Liao  (1990),  will  later  be  verified  in  a  particular  example. 


Average  micropolar  stress  and first  moment  of  stress 


Because  Eq.  43  applies  to  arbitrary  values  of  fa  and  cijk,  the  following  relations  are 
obtained: 


fai  +eiklfalk  =^{eikJic{xbM -4)~Xj(xck  -xak))+m'(xbj -xaj))^^Mfx' 

V  CG 1  *  eeE 


ij- 1,2,3 

^kji +  faki  =  77 fi  (xjXk  -XjXk)  =—£fi  Xj  xk  i ,j,  k - 1 ,2,3 

V  cel  * 


(2.49) 

(2.50) 


eeE 


Therefore  the  external  moments  M°c ,  which  result  from  external  contact  force  fe  and/or 
contact  moment  m‘ ,  generate  not  only  asymmetric  stress  but  also  couple  stress  and  first 
stress  moment.  This  result  is  in  agreement  with  Eq.  26,  which  states  that  couple  stresses 
are  required  to  balance  asymmetric  stresses.  However,  the  present  approach  provides 
only  the  sums  /I),  +  e,fLjlk  and  T.kJi  +  T,Jki ,  and  unfortunately  not  each  term  ]Jj,  and  . 


Alternate  definition  of  average  stress 

The  average  stress  can  also  be  defined  based  on  statics,  instead  of  virtual  work  (e.g., 
Cundall  and  Strack,  1978).  The  average  stress  within  volume  Fis  defined  as  the 
weighted  average  of  the  stress  o°  for  each  particle  a  of  B : 


'  oeB 

where  Va  is  the  volume  of  particle  a,  and: 


(2.51) 


(2.52) 
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By  using  Eq.  34,  and  replacing  the  traction  vector  Tt  with  discrete  contact  force  f-  ,  o° 
becomes: 


(2.53) 


Because  the  contact  forces  are  opposite  at  internal  contacts  (i.  e.,  fja  =  -ff ),  the 
average  stress  if*  is: 


v  azBceC.  "  cel  V  ezE  Y  ezE 


Note  that  x‘  in  Eq.  54  refers  to  contact  point  e,  whereas  x°e  in  Eq.  44  refers  to  the  center 
of  the  particle  a  where  contact  e  takes  place.  cfy  and  a y  are  related  through: 


(2.55) 


The  symmetry  of  ojj  results  from  the  equilibrium  of  moments  about  the  coordinate 
origins  for  particle  a ,  i.e.: 


^ijk® jk 


ceCa 


(2.56) 


The  stress  o  *  is  therefore  symmetric  because  it  is  the  weighted  sum  of  symmetric  <7° . 
The  symmetry  of  if*  can  also  be  shown  by  using  Eq.  54  and  invoking  the  equilibrium  of 
external  moments  about  the  coordinate  origin  (i.e.,  Eq.  22). 


The  symmetry  of  cr  *  has  significant  implications  in  computational  granular  mechanics, 

especially  for  the  computer  simulations  using  dynamic  relaxation  to  solve  the  equilibrium 
equations  of  statics  (e.g.,  Cundall  and  Strack,  1978;  Bardet  and  Proubet,  1991).  When 
there  is  no  moment  at  contacts,  a *  should  be  symmetric,  and  any  computed  asymmetry 

of  o'  should  be  interpreted  as  inaccurate  calculation  and/or  lack  of  static  equilibrium. 
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In  the  case  of  spherical  and  cylindrical  particles  of  radius  Ra,  xct  =  x°  +  Ran ■ ,  the  average 
stress  o'j  becomes: 

t-y'Ly  iw  +  4X*-  X  (2'57) 


aeBv  a  ceCa 


V  “  V  n 

v  aeB  v a  ceCa 


aeB  ceCa 


Equation  57  is  the  same  as  that  obtained  by  Cundall  and  Strack  (1979). 

Examples 

We  will  illustrate  the  circumstances  of  stress  asymmetry  in  the  case  of  double  and 
multiple  layer  interfaces  filled  with  cylindrical  or  spherical  particles. 

Example  1:  Double  layer  interface 

The  stresses  atj  and  a*  can  be  calculated  analytically  for  the  particular  example  of  Fig. 

3,  which  represents  an  interface  made  of p  columns  each  having  two  particles.  The 
columns  have  the  same  height  but  are  made  of  particles  of  various  diameters.  The  particle 
assembly  is  subjected  to  the  normal  force  Np  and  shear  force  Sp,  which  are  assumed  to  be 
distributed  evenly  onto  each  particle  column;  the  normal  and  shear  forces  acting  at  all 
contacts  are  N=  N/p  and  S  =  S/p,  respectively.  The  contact  direction  between  two 
particles  is  identical  for  all  particles  in  contact.  It  is  characterized  by  the  unit  vector  n  of 
component «/  =  sin  G  and  =  cos0.  The  equilibrium  of  forces  and  moments  for  all 

particles  and  the  top  plate  are  satisfied  when: 


S  =  Nnf{l  +  n2)=  N  tan  (6/2) 


(2.58) 


All  the  contact  forces  have  the  same  inclination  6/2  relative  to  the  contact  direction  n. 
The  contacts  do  not  slip  and  the  interface  remains  stable  as  long  as  6  remains  smaller 
than  2<j),  where  <j)  is  the  friction  angle  between  two  particles.  For  the  calculation  of  <7y,  it 

is  convenient  to  select  the  coordinate  axis  at  the  center  of  particle  2.  In  this  coordinate 
system,  the  center  coordinates  of  particle  1  are: 


x,le  =  nf  and  x\e  =  n2L 


(2.59) 
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where  L  =  Ri  +  R2.  By  using  Eqs.  44  and  59,  the  stresses  <7tJ  are: 


°n=-"iS/A.  on—$*S&SlA.  (260) 

<J12  =  -(l  +  n2  )S/ A ,  and  <t2]  =  -n2  S/ A 

where  A  =  V/L,  and  V  is  the  average  volume  of  particle  columns.  For  interfaces  filled 
with  cylindrical  or  spherical  particles,  V can  be  evaluated  as  follows: 

V  =  {\  +  n2)LWxW2  (2.61) 

where  Wi  is  the  average  width  of  particle  columns  in  the  ^/-direction.  For  cylindrical 
particles,  W2  is  the  particle  length.  For  spherical  particles,  W2  is  the  average  width  of 
particle  columns  in  the  ^-direction.  Both  fVj  and  W2  depend  on  the  density  of  particles  in 
the  interface. 

The  external  moments  M°e  acting  on  particles  1  and  2  are: 


M\  =  RXS  and  M\  =  R2S 

By  using  Eq.  47  or  60,  the  amplitude  of  stress  asymmetry  is: 


(2.62) 


G\1  °21 


.MhMl  =-s/a 


(2.63) 


The  stress  asymmetry  vanishes  when  S  =  0,  i.e.,  when  the  particle  columns  are  vertical. 
As  previously  stated  in  Eq.  48,  Oy  is  not  necessarily  symmetric  for  particles  of  identical 

radius  (e.g.,  Ri  =  R2).  The  couple  stress  and  first  stress  moment  is: 


/i13  +  Z,2|  —  E|,2  —  RxH\  S/ A  and  fi2 3  +  "Z22]  E2]2  —  /?,m2  S/A  (2.64) 


For  the  calculation  of  <7* ,  it  is  convenient  to  select  the  coordinate  axis  at  the  lowest 
external  contact.  The  coordinates  of  the  highest  contact  are  therefore: 


jcf  =  nxL  and  x2  =  (1  +  n2)L 
and  the  stress  <7*  is: 


(2.65) 


c t,*,  =  5'/y4,  <722  =  A,  and <rl2  =  cr21  =  -(l  +  n2)S/A 


(2.66) 


As  expected,  <7,*  is  symmetric.  a ‘and  <7,-,  are  related  through: 

axx  =  an,  cx 22  =  a22  -  ^  +  n^S/At  c'n  =  a2X,  and  o2X  =  ct2I  - 5/^  (2.67) 

n\ 

Figure  4  compares  the  normalized  variation  of  tf’and  a tj  when  the  angle  6  =  tan''(«//«^) 
varies  from  0  to  90°.  The  stress  asymmetry  (i.e.,  on  -  a2X )  remains  constant  and 
independent  from  6. 

Example  2:  Multi-layer  interface 

As  shown  in  Fig.  5,  the  multi-layer  interface  model  is  made  of  p  columns  of  particles, 
each  column  /  having  qt  particles  (q,  must  be  an  even  number).  This  multi-layer  model  is 
a  generalization  of  the  two-layer  model,  which  corresponds  to  q,  =  2  for  i  =  1  ,...,p.  The 
particle  columns  have  the  same  height  and  contact  direction  n,  but  may  have  different 
numbers  qt  of  particles.  All  the  particles  and  the  top  platen  are  in  static  equilibrium  when 
Eq.  58  is  satisfied;  the  contact  forces  are  identical  to  those  of  the  two-layer  model.  The 
average  stresses  in  the  interface  can  be  determined  by  averaging  the  average  stresses  in 
each  column.  Hereafter,  we  will  calculate  only  the  average  stresses  in  a  column. 

After  selecting  the  coordinate  axis  at  the  center  of  particle  qt,  the  center  coordinates  of 
particle  1  are: 


x,  =  nxL  and  x2  =  (1  +  n2)L  -  (Rx  +  Rq  ) 


(2.68) 


where  L  =  . .  The  average  stresses  cy  in  column  i  are: 

7=1 


^11  —  ®22  ~ 


(l  +  n2f  ,  *1  +  Rq,  1  +  »2 


-S/A  +  - 


R  +  R 

<JX2  =  -(l  +  n2)S/A,  and  a2,  =  -(l  +  n2)S/A  +  '  -  S 


V  n, 
R, 

+-  — 

V 


S , 


(2.69) 


-  Page  54  - 


where  A  and  V  are  defined  as  for  the  double  layer  interface.  The  values  of  <7U  and  an 
are  identical  to  those  in  Eqs.  60,  whereas  the  corresponding  values  of  o22  and  cf21  are 
different.  The  differences  however  vanish  when  The  stresses  <7 *  in  the  multi¬ 

layer  are  identical  to  those  of  the  double  layer  interface  (i.e.,  Eq.  66).  As  predicted  by  Eq. 
69,  Oy  converges  toward  (7*  when  volume  V becomes  large  (i.e.,  q  — >  °°).  The  external 

moments  M*  acting  on  particle  1  and  q,  are: 


M\  =  RXS  and  Mf  =  Rq  S 

These  external  moments  are  responsible  for  the  stress  asymmetry  as  follows: 


(2.70) 


an-o2l=-S(Rl+Rqi)/V 


(2.71) 


The  amplitude  of  stress  asymmetry  decreases  with  the  number  of  particles  in  the 
columns,  and  increases  with  the  size  of  particles  at  the  top  and  bottom  of  columns.  The 
stress  becomes  symmetric  when  the  column  height  becomes  infinite  (i.e.,  q  — >  °°).  The 
couple  stress  and  first  stress  moment  are: 


/^I3  +  ^121  ^112  -  R\n\  3/ -A 

+  t22i  - E2I2  =  i?,(l  +  n2)S/A - Rl(R]  +  RJS/V 


(2.72) 


In  contrast  to  the  asymmetric  stress  components,  fiy  +  does  not  decrease  when  the 

volume  V becomes  large  (i.e.,  q  — >  °o). 
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Figure  21. 


Multi-layer  interface  model  for  calculation  of  average  stress. 
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Discussion 


The  asymmetry  of  stress  depends  on  the  way  stresses  are  defined.  Stresses  defined  by 
statics  are  symmetric  when  there  is  no  contact  moment.  However  stresses  defined  by 
virtual  work  may  become  asymmetric  when  there  is  no  contact  moment.  Our  analysis 
differs  from  previous  ones  (e.g.,  Christoffersen  et  al.,  1981;  and  Chang  and  Liao,  1990) 
because  we  considered  external  moments,  and  established  the  circumstances  and 
amplitude  of  stress  asymmetry.  In  agreement  with  Caillerie  (1991),  we  found  that  the 
asymmetry  originates  from  external  moments,  and  that  the  amplitude  of  stress  asymmetry 
decreases  with  the  size  of  the  granular  volume.  The  stress  asymmetry  is  therefore  more 
detectable  in  elongated  samples  subjected  to  external  moments  on  their  boundary.  Bulky 
samples  subjected  to  small  external  moments  are  likely  to  display  negligible  stress 
asymmetry.  The  stress  asymmetry  can  rightfully  be  neglected  in  large  masses  of  granular 
media  far  away  from  boundaries  with  external  moments.  However,  it  may  become 
important  in  interfaces  with  significant  external  moments.  There  is  a  need  for  verifying 
these  findings  through  computer  simulations  and  laboratory  experiments. 

Conclusion 

We  have  derived  the  conditions  for  the  asymmetry  of  stress  in  granular  materials,  and 
shown  that  there  is  asymmetry  even  when  the  particle  contacts  do  not  transmit  moments. 
This  stress  asymmetry  is  obtained  when  the  stress  is  defined  from  virtual  work,  but  is  lost 
when  the  stress  is  defined  from  statics.  The  asymmetry  results  from  external  moments 
applied  by  external  forces.  We  also  show  that  the  amplitude  of  stress  asymmetry 
decreases  with  the  ratio  V/S  between  surface  S  and  volume  V.  When  V/S  becomes  very 
large,  the  stress  asymmetry  disappears,  with  or  without  external  contact  moments. 
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PART  III.  EXPERIMENTAL  INVESTIGATION 

Parts  I  and  II  of  this  report  covered  some  of  the  computational  and  theoretical  approaches 
for  investigating  material  instability  in  granular  materials.  Part  III  devises  an 
experimental  approach  to  provide  these  computational  and  theoretical  approaches  with 
relevant  and  detailed  experimental  data  sets.  Specimens  of  idealized  granular  media  were 
constructed  and  tested  in  the  laboratory,  and  the  displacement  and  rotation  of  particles 
were  measured  using  stereophotogrammetry.  The  laboratory  tests  were  conducted  at  the 
University  of  Southern  California;  the  stereophotogrammetric  measurements  were  carried 
out  at  the  Joseph  Fourier  University,  in  Grenoble,  France. 

Stereophotogrammetry  was  found  to  be  the  best  suited  techniques  for  measuring 
accurately  the  displacements  and  rotations  of  a  large  number  of  nearly  identical  particles. 
At  the  time  of  this  investigation,  this  optical  technique  revealed  to  be  more  accurate  and 
practical  than  computer  vision  methods  based  on  existing  digital  cameras.  The  accurate 
experimental  results  generated  in  this  study  are  useful  to  assess  the  applicability  of 
higher-order  continua  to  granular  behavior,  and  are  therefore  instrumental  to  understand 
material  instabilities  within  granular  media. 

The  first  section  of  Part  III  reviews  the  background  that  motivated  our  experimental 
investigation.  The  second  section  describes  the  experimental  set-up  and  the 
stereophotogrammetric  measurements.  The  third  and  last  section  describes  experimental 
data  processing  and  presents  processed  results. 

Background 

Deformations  within  granular  soils  are  commonly  concentrated  in  narrow  zones  called 
shear  bands,  the  thickness  of  which  are  8-10  times  the  mean  grain  diameter  (Roscoe 
1970;  Muhlhaus  and  Vardoulakis  1987).  The  relation  between  shear  band  width  and 
grain  size  has  profound  implications  in  laboratory  and  centrifuge  testing  of  reduced-scale 
models.  It  implies  that  grain  sizes  must  also  be  scaled  down  in  laboratory  and  centrifuge 
models  for  simulating  the  progressive  failure  of  in-situ  granular  masses. 
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Hill  (1962),  Mandel  (1963),  Rudnicki  and  Rice  (1976)  successfully  analyzed  the 
emergence  and  inclination  of  shear  bands,  but,  since  their  constitutive  theories  included 
no  internal  length  scale,  they  did  not  address  the  problem  of  shear-band  thickness. 

Aifantis  (1984, 1987)  and  Zbib  and  Aifantis  (1989)  introduced  an  internal  length  scale  by 
including  second-  and  fourth-order  strain  gradients  into  the  yield  condition  of  plasticity. 
The  higher-order  strain  gradients  allowed  them  to  account  for  the  thickness,  spacing,  and 
velocity  of  shear  bands  in  metals.  Vardoulakis  and  Aifantis  (1989)  applied  a  similar 
approach  to  soils  by  adding  second-order  strain  gradients  to  the  flow  rule  of  plasticity. 

As  an  alternative  to  higher-order  strain  gradient  theories,  Muhlhaus  (1986),  Muhlhaus 
and  Vardoulakis  (1987),  and  Vardoulakis  (1989)  introduced  an  internal  length  scale  by 
using  the  micropolar  (Cosserat)  theory.  Although  based  on  different  physical 
assumptions,  both  micropolar  and  strain  gradient  approaches  successfully  detect  the 
emergence  of  shear  bands  and  calculate  their  inclination  and  thickness.  Zbib  and  Aifantis 
(1989)  also  analyzed  the  evolution  of  the  deformation  within  the  shear  bands  of  metals. 
They  could  not  apply  it  to  soils  due  to  lack  of  relevant  experiments  on  soils. 

The  structure  of  shear  bands  in  real  granular  materials  is  difficult  to  investigate  by  using 
laboratory  experiments.  The  radiographic  techniques  with  lead  shots  (Roscoe,  1963) 
estimate  the  shear  band  thickness  but  are  not  accurate  enough  to  measure  the  localized 
field  of  displacement.  The  stereophotogrammetric  technique  of  Desrues  (1984)  measures 
the  displacement  but  not  the  rotation  of  particles.  The  numerical  simulation  techniques  of 
Bardet  and  Proubet  (1989,  1991,  1992),  and  Bardet  (1994)  provided  extensive 
information  on  the  displacement  and  rotation  of  particles  but  are  only  numerical 
substitutes  to  actual  laboratory  experiments. 

Herein,  we  present  an  experimental  approach  that  generates  experimental  data  sets  on 
cylindrical  rods  and  is  useful  to  investigate  the  application  of  higher-order  continuum 
theories  to  granular  media.  Cylindrical  rods,  in  place  of  real  granular  material,  are  tested 
under  confined  axial  compression,  in  an  attempt  to  emulate  plane-strain  compression  tests 
on  real  soils.  Stereophotogrammetry  (Desrues,  1984)  are  used  to  determine  the 
displacement  and  rotation  of  particles.  Stereophotogrammetry  was  preferred  over  other 
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techniques  for  measuring  the  displacement  of  particles.  It  visualizes  displacement  fields 
as  a  false  relief,  and  automatically  tracks  particles  between  two  successive  photographs. 
This  is  a  great  benefit  in  the  case  of  large  assemblies  of  nearly  identical  particles. 
Stereophotogrammetry  automatically  digitizes  and  recognizes  individual  particles 
between  successive  photographs. 

Axial  compression  of  idealized  granular  media 

Sample  composition  and fabrication 

As  listed  in  Tables  1  and  2  and  shown  in  Fig.  23,  two  tests,  hereafter  referred  to  as  A  and 
B,  were  carried  out;  the  specimens  of  tests  A  and  B  were  made  of  a  total  of  1848 
cylindrical  acrylic  rods  having  for  nominal  diameter  4,  6  and  8  mm.  Specimens  A  and  B 
had  a  slightly  different  slenderness  ratio.  All  rods  are  made  of  transparent  acrylic  material 
with  an  average  density  of  1 .3  g/cm3.  The  rods  were  carefully  cut  to  the  same  length  of 
10  cm;  their  faces  were  perpendicular  and  transparent  as  much  as  possible.  Front  faces 
were  half  painted  black  in  order  to  track  simultaneously  rotations  and  displacements. 
After  a  careful  inspection,  the  acrylic  rods  were  found  to  be  slightly  imperfect  cylinders. 
Their  diameters,  as  measured  with  a  caliper,  slightly  vary  in  the  radial  and  longitudinal 
directions.  Table  3  summarizes  average  diameter  and  diameter  variation  A<j>  for  a  small 
sample  of  arbitrarily  selected  rods.  The  average  diameter  is  8.12,  6.19  and  4.22  mm,  in 
comparison  to  the  nominal  values  of  8  mm,  6  mm  and  4  mm.  Measured  diameters  vary 
typically  less  than  1%.  The  black  paint  was  also  found  not  to  divide  exactly  particles  in 
half,  and  could  not  be  used  to  determine  particle  diameters. 

Table  1.  Particle  size  distribution  for  samples  A  and  B  of  idealized  granular  material. 


Nominal 

Diameter 

(mm) 

Measured 

Average 

Diameter 

(mm) 

Number  of  Particles 

Complete 

Specimen 

Stereo  Sample 
A 

Stereo  Sample 

B 

4.0 

4.22 

1040 

458 

419 

6.0 

6.19 

532 

243 

248 

8.0 

8.12 

266 

129 

120 

-  Page  61  - 


* 

Figure  23.  Assembly  of  cylindrical  rods  subjected  to  axial  compression. 


Table  2.  Result  summary  of  tests  A  and  B 


Slenderness 

Initial  Young’s 

Ultimate 

Friction 

Strain  at 

Number 

TEST 

ratio 

Modulus 

Strength 

angle  at 
peak  failure 

peak 

failure 

of 

particles 

H/L 

E  (MPa) 

C^max 

(kPa) 

^max  (deg) 

£$max 

(%) 

digitized 

A 

2.43 

304 

84 

18.4 

B 

2.81 

574 

88 

18.3 

3.24 

787 

Table  3.  Variation  of  rod  diameters  on  random  samples. 


Diameter  ^mm) 


Trial  No. 

1 

2 

3 

4 

5 

Average 

Min 

Max 

A«k  (mm)  A4>/<t>  (%) 

Average 

8  mm  rod  1 

8.26 

8.18 

8.19 

8.23 

8.22 

8.18 

8.26 

0.08 

0.97 

8.12 

2 

8.07 

8.11 

8.15 

8.09 

8.05 

8.09 

8.05 

8.15 

0.10 

1.24 

3 

8.04 

8.05 

8.08 

8.07 

8.02 

8.05 

8.02 

8.08 

0.06 

0.75 

6  mm  rod  1 

6.11 

6.08 

6.08 

6.12 

6.07 

6.09 

6.07 

6.12 

0.05 

0.82 

6.19 

2 

6.31 

6.26 

6.27 

6.30 

6.26 

6.28 

6.26 

6.31 

0.05 

0.80 

3 

6.20 

6.21 

6.15 

6.20 

6.19 

6.19 

6.15 

6.21 

0.06 

0.97 

4  mm  rod  1 

4.25 

4.25 

4.22 

4.21 

4.24 

4.23 

4.21 

4.25 

0.94 

4.22 

2 

4.24 

4.23 

4.22 

4.26 

4.25 

4.24 

4.22 

4.26 

0.94 

3 

4.19 

4.17 

4.22 

4.19 

4.21 

4.20 

4.17 

4.22 

0.05 

«  A  A  ,  A  /O/.  \  - 

1.19 

ft  nc 

AverageA<(>/<l)(%)  =  0.96 


4 
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Experimental  setup  for  axial  compression 

Figures  24  shows  the  experimental  set  up.  The  test  specimen  of  Fig.  23  is  enclosed  in  a 
0.15  mm  thick  transparent  latex  membrane  and  has  circular  loading  platens  at  its  top  and 
bottom.  The  membrane  is  clamped  to  both  platens  with  compression  rings.  During  the 
sample  fabrication,  the  membrane  is  stretched  on  a  rectangular  mold,  and  the  rods  are 
manually  positioned  inside  the  mold  with  their  axis  parallel  to  one  another.  Specimens  A 
and  B  have  therefore  a  different  packing,  in  addition  to  a  different  slenderness  ratio. 
Vacuum  is  applied  inside  the  membrane,  which  is  equivalent  to  applying  an  external 
constant  confining  pressure  to  specimens.  The  equivalent  confining  pressure  was  equal  to 
94.8  kPa.  The  axial  compression  is  applied  by  raising  the  bottom  platen  at  the  constant 
rate  of  5.6  mm/min. 

As  shown  in  Fig.  24,  a  light  box  behind  the  sample  reflects  and  diffuses  the  light  of  a  120 
watts  bulb  through  a  polished  glass.  A  Nikon  FE2  35  mm  camera  is  positioned 
approximately  one  meter  in  front  of  the  specimen.  It  is  equipped  with  a  macrolens  to 
limit  optical  distortion.  The  film  was  a  hard  contrast,  black  and  white  ortho  film,  25  ASA, 
with  a  fine  grain.  The  camera  aperture  speed  was  set  to  one  second. 


Reaction 

frame 


Polished 


'////fSS//S/S//S/SS/SSSSSSSSSSSSS*SSSSSSSSSSSSSSS*fSSSS*SSSfSSSSfSS'+*sffsrsslrs*r*srf  . .  ' 

Figure  24.  Experimental  setup  for  axial  compression  of  idealized  granular  materials. 
Results 


Specimens  A  and  B  were  tested  similarly  by  increasing  the  axial  force  while  keeping  the 
confining  pressure  (i.e.,  internal  vacuum)  constant.  The  stress-strain  curves  obtained 
during  tests  A  and  B  are  shown  in  Fig.  25.  The  axial  strain  e  is 


£  =  Ah  /  h0  (3.1) 

where  Ah  is  the  vertical  displacement  of  the  lower  platen,  and  h0  is  the  initial  sample 
height.  The  axial  stress  c  is 

o=F/A0  (3.2) 

where  F  is  the  measured  axial  load,  and  A  0  is  initial  cross-sectional  area.  Figure  25 
shows  the  stress-train  response  until  6%  axial  strain.  Both  samples  failed  at 
approximately  at  3%  strain.  Table  2  summarizes  the  calculated  results  for  initial  Young 
modulus,  peak  and  residual  friction  angles,  and  strain  at  peak  failure  for  tests  A  and  B. 
The  symbols  in  Fig.  25  indicate  the  axial  strains  at  which  the  photographs  were  taken. 


Eight  photographs  were  digitized  in  tests  A  and  B.  Four  photographs  were  taken  at  points 
A1  to  A4  during  test  A,  and  at  points  B 1  to  B4  during  test  B. 


♦ 


0 


0 


0  2  4  6 

Strain  (%) 

Figure  25.  Axial  stress-strain  response  of  samples  A  and  B,  and  location  of  photographs. 
Stereophotogrammetry  Measurement 

The  eight  photographs  taken  during  tests  A  and  B  were  processed  by 
stereophotogrammetry,  an  optical  technique  especially  suited  for  tracking  the  motion  of  a 
large  number  of  nearly  identical  particles. 

Principle  of  stereophotogrammetry  for  measurement  of  displacement 

The  principles  of  stereophotogrammetric  measurement  of  displacement  fields  are 
described  in  Desrues  (1984).  The  stereophotogrammetry  equipment  is  called  stereograph. 
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The  absolute  error  on  measured  displacement  varies  with  quality  of  negatives  and 
operator  experiences.  In  optimum  conditions,  the  stereograph  can  measure  1  pm  at  the 
negative  scale,  i.e.,  as  small  as  the  grains  coating  negatives.  In  our  experiment,  the 
absolute  error  is  assumed  to  be  in  the  order  of  10  pm  as  we  used  35  mm  negatives, 
instead  of  larger  negatives. 

Application  of stereophotogrammetry  to  kinematics  of  particles 

The  35  mm  negatives  taken  during  tests  A  and  B  were  enlarged  to  obtain  12.7  x  17.8  cm 
(5x7  inch)  negatives,  the  stereograph  using  only  transparent  media.  The  negatives 
corresponding  to  two  successive  photographs  (i.e.,  A1  and  A2)  were  placed  on  the 
stereograph  platens.  At  first,  three  fixed  reference  points  (FI,  F2  and  F3  in  Fig.  26)  were 
selected  close  but  outside  the  sample.  These  three  points,  which  correspond  to  fixed 
objects  in  the  laboratory,  are  references  for  defining  a  common  coordinate  system  to  all 
photographs.  In  principle,  the  stereograph  should  indicate  that  fixed  points  do  not  move 
(i.e.,  move  only  within  the  order  of  machine  accuracy).  For  experimental  reliability,  these 
points  are  digitized  twice.  Second,  the  stereograph  software  asks  to  define  the  area  under 
investigation  by  selecting  four  points  (see  Fig.  26). 

Positioning  particles  requires  that  three  points  be  digitized  for  yielding  center  position 
and  radius  of  individual  particle.  Three  points  are  needed  because  the  black  paint  does 
not  exactly  delineate  particle  diameters.  As  shown  in  Fig.  27,  points  A  and  B  are  obvious 
candidates  as  they  are  easily  identified.  However,  point  C  is  a  less  obvious  choice.  As 
explained  in  a  later  section  on  error  minimization,  point  C  is  chosen  as  the  point  farthest 
away  from  points  A  and  B. 


ABSTRACT 


The  research  aims  at  investigating  the  physical  origins  of  shear  band  instability  in 
particulate  media,  and  the  applicability  of  higher-order  continuum  theories  for  describing 
material  instability  in  granular  media.  The  particular  research  objectives  are  (1)  to  review 
the  computational  tools  used  for  simulating  shear  bands  in  idealized  granular  media,  (2) 
to  explore  a  critical  micro-macro  transition  relevant  to  material  instability  (i.e.,  stress 
symmetry),  and  (3)  to  devise  laboratory  experiments  for  generating  some  experimental 
data  sets  useful  for  numerical  and  theoretical  investigations  of  material  instability  in 
granular  media.  Our  methodology  combines  computational  granular  mechanics,  higher- 
order  continuum  mechanics,  and  laboratory  experiments. 

Our  comprehensive  review  of  past  work  on  computational  granular  mechanics  reveals  the 
diversity  of  numerical  codes  in  granular  mechanics,  and  the  large  number  of  assumptions 
entering  these  codes. 

Our  theoretical  investigation  shows  that  the  average  stress  in  granular  media,  as  defined 
from  virtual  work,  may  be  asymmetric  in  the  absence  of  contact  moments.  We  specify  the 
circumstances  and  amplitude  of  stress  asymmetry,  and  calculate  the  corresponding  couple 
stress  and  first  stress  moment.  We  also  show  that  the  average  stress  is  always  symmetric 
when  it  is  alternately  defined  by  using  statics  and  no  contact  moment.  The  stress 
asymmetry,  which  results  from  external  moments,  decreases  with  the  volume  size.  The 
asymmetric  stress,  couple  stress  and  first  stress  moment  are  analytically  calculated  in 
particular  examples  relevant  to  shear  band  instability. 

Our  experimental  methodology  is  based  on  the  application  of  stereophotogrammetry  to 
the  kinematics  of  assembly  of  cylindrical  rods,  simulating  granular  media.  We  are  able  to 
measure  accurately  the  kinematics  of  a  large  number  of  nearly  identical  particles  in  the 
laboratory.  Our  experimental  methodology  has  generated  new  experimental  data  sets  on 
granular  media  useful  to  re-examine  the  applicability  of  higher-order  continua  to  granular 
media,  and  are  therefore  instrumental  to  understand  material  instability  in  granular  media. 


100  120  140  160  180  200 


(X) 

Figure  26.  Geometry  of  sample  A,  the  control  points  FI  F2  and  F3  and  the  digitizing 
area. 
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Figure  27.  Calculation  of  center  position  and  radius  of  particle  using  three  different 
methods:  (a)  Method  A,  (b)  Method  B,  and  (c)  Method  C. 

Determination  of  Assembly  Geometry 

The  stereograph  measures  the  displacement  of  points  A,  B  and  C  for  individual  particles 
between  two  successive  photographs.  As  shown  in  Fig.  27,  points  A  and  B  are  about 
diametrically  opposite,  and  at  the  boundary  of  the  painted  area.  Points  A  and  B  are  used 
to  track  particles  between  successive  photographs.  Point  C  is  used  to  define  the  center 
and  radius  of  particles,  assuming  that  particles  are  perfectly  rigid  and  circular.  The 
stereograph  cannot  track  point  C  because  there  is  not  special  mark  on  particles  associated 
with  point  C. 

Assumption  of  rigid  particles 

The  assumption  of  rigid  particles  was  tested  by  measuring  the  variation  AAB  of  distance 
AB  between  the  successive  steps  A1-A2,  A2-A3  and  A3-A4.  In  theory,  AAB  should  be 
equal  to  zero  if  particles  were  purely  rigid  and  measurements  perfectly  accurate.  As 
shown  in  Fig.  28,  the  distribution  of  MB  indicates  that  AB  varies  less  than  1  %  between 
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two  successive  steps. 
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Figure  28.  Variation  of  distance  AB  at  different  axial  compression  of  sample  A. 
Particle  center  and  radius 

As  shown  in  Fig.  27,  three  different  methods,  labeled  A,  B  and  C,  were  used  to  calculate 
the  particle  center  and  radius  from  the  raw  data  measured  by  the  stereograph.  The 
corresponding  calculation  steps  are  detailed  in  the  FORTRAN  programs  GET  AB  in 
APPENDIX. 


Method  A 

Method  A  calculates  the  center  position  and  radius  of  particles,  assuming  that  points  A,  B 
and  C  are  on  a  circle.  As  shown  in  Fig.  27a,  the  circle  center  O  is  the  intersection  of  lines 
PQ  and  RS  passing  through  the  center  of  segments  AB  and  AC  and  perpendicular  to  AB 
and  AC,  respectively.  Point  O,  which  is  the  intersection  of  lines  PQ  and  RS,  has  the 
following  coordinates  (xo,yo): 


(x^  +  (yA  -TbVo  —(?a  ~xb +(?a 

- 

\(xA -XC)xo  +G ’a  “Tc)fo  =  (*.4  ~xc)l'^+(^A  -yc)!^ 


(3.3) 
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where  (xA,yA),  (xb,  Yb)  and  (xc,yc)  are  the  coordinates  of  points  A,  B,  and  C, 
respectively.  Coordinates  (x0,yo)  are  always  uniquely  defined  provided  that  points  A,  B 
and  C  are  not  aligned.  The  particle  radius  R  is  the  distance  OA,  which  is  the  same  as  OB 
and  OC. 

Methods  B  and  C 

As  shown  in  Fig.  27b,  Method  B  assumes  that  distance  OD  between  segment  AB  and 
point  O  remains  constant  for  a  given  particle.  Method  B  uses  points  A  and  B;  it  only  uses 
point  C  at  the  first  step  to  calculate  the  initial  value  of  OD. 

As  shown  in  Fig.  27c,  Method  C  assumes  that  angle  OAB  remains  constant  for  a  given 
particle.  Like  Method  B,  Method  C  uses  only  point  C  at  the  first  step  to  calculate  the 
initial  value  of  angle  OAB.  Like  methods  A  and  B,  method  C  does  not  fix  the  distance 
AB.  Methods  A,  B  and  C  were  found  to  give  identical  results. 

Results  on  particle  diameter 

In  theory,  the  measured  diameters  of  rigid  particles  should  remain  constant  between 
successive  photographs.  In  practice,  the  particle  diameters,  as  calculated  by  methods  A,  B 
and  C,  were  found  to  vary  slightly  between  steps.  Figures  29  and  30  show  the  cumulative 
distributions  of  rod  diameter,  which  were  calculated  from  the  stereograph  at  the  four 
stages  A1  to  A4,  and  B1  to  B4.  As  shown  for  sample  A  in  Figs.  29  and  30,  the  four 
curves  labeled  A1  to  A4,  and  B1  to  B4,  are  in  good  agreement,  which  establishes  that 
Method  A,  B  and  C  gives  very  similar  rod  diameter  between  successive  photographs. 

The  staircase  line  in  Figs.  29  and  30,  which  is  labeled  “total”,  represents  the  cumulative 
distributions  of  rod  diameter,  assuming  only  three  different  diameters.  The  smooth 
measured  distributions  indicate  that  particle  diameters  vary  according  to  a  tri-modal 
statistical  distribution  centered  about  the  average  values  of  Table  3.  Figures  29  and  30 
also  show  that  samples  A  and  B  and  the  complete  specimen  have  slightly  different 
distribution  of  particles.  Not  all  particles  of  specimens  A  and  B  were  digitized.  Particles 
close  to  the  loading  caps  and  close  to  the  lateral  boundaries  were  excluded.  As  indicated 
in  Table  1,  830  particles  were  digitized  in  Specimen  A,  and  787  particles  in  Specimen  B. 


Percent  Finer 


As  shown  in  Table  1,  the  corresponding  number  of  particles  with  different  diameter  in 
Samples  A  and  B  were  calculated  from  Figs.  29  and  30. 


Figure  29.  Cumulative  rod  diameter  distribution  curves  for  sample  A  based  (1)  on 

stereo  measurement  at  axial  compression  stages  Al,  A2,  A3,  and  A4,  and  (2) 
on  the  complete  specimen  composition  in  Table  1 . 


Figure  30.  Cumulative  rod  diameter  distribution  curves  for  sample  B  based  (1)  on 

stereo  measurement  at  axial  compression  stages  Bl,  B2,  B3,  and  B4,  and  (2) 
on  the  complete  specimen  composition  in  Table  1. 

Figure  3 1  shows  that  there  is  little  difference  between  the  cumulative  distribution  curves 
for  the  complete  specimen  and  Samples  A  and  B.  Figure  32  shows  the  variation  of 
particle  diameters  which  is  calculated  between  the  axial  compression  stages  A2,  A3  and 
A4  relatively  to  stage  A 1.  Figure  31  implies  that  measured  particle  diameters  vary  with 
axial  strain.  These  small  variations  result  from  a  combination  of  experimental  errors  and 
assumptions  on  particle  rigidity  and  shape. 


Figure  3 1 .  Cumulative  rod  diameter  distribution  curves  for  the  complete  specimen  and 
Samples  A  and  B. 


Figure  32.  Variation  of  particle  diameter  as  measured  by  stereograph  at  axial 
compression  stages  A2,  A3,  and  A4  relatively  to  stage  Al. 


Error  on  center  position  and  radius 

Figure  33  shows  the  graphical  determination  of  errors  associated  with  the  determination 
of  particle  center  O.  This  graphical  construction  is  more  convenient  than  algebraic 
derivations  of  error  based  on  Eq.  3.3.  The  possible  positions  of  point  O  is  encompassed 
by  the  shaded  area,  which  is  generated  when  points  A,  B  and  C  move  independently 
within  the  experimental  error  Ax  (i.e.,  when  points  A,  B  and  C  describe  the  circles  of 
radius  Ax  centered  about  points  A,  B  and  C).  This  graphical  construction  proves  that  the 
error  depends  of  the  relative  position  of  points  A,  B  and  C.  When  points  A  and  B  are 
approximately  diametrically  opposite,  the  smallest  error  is  obtained  when  C  is  the 
furthest  away  of  points  A  and  B  (see  Fig.  33a).  When  point  C  gets  close  to  point  A  or  B, 
the  shaded  area  becomes  elongated,  and  the  error  increases  significantly  in  the  direction 
perpendicular  to  AB  (see  Fig.  33b).  In  summary,  an  accurate  determination  of  particle 
center  and  radius  requires  to  take  point  C  as  far  as  possible  from  points  A  and  B.  In  the 
present  investigation,  the  optimum  condition  was  selected;  the  error  associated  with  the 
position  of  point  O  is  comparable  to  those  of  points  A  and  B. 


Figure  33.  Definition  of  error  for  particle  center:  (a)  Point  C  far  away  from  points  A  and 
B,  and  (b)  Point  C  close  to  points  A  or  B. 
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Particle  angular  position 


» 


As  sketched  in  Fig.  34a,  when  the  particle  remains  rigid,  the  particle  rotation  about  its 
center  is  identical  to  the  rotation  of  segment  AB.  Thus,  the  angular  position  0  of  the 
particle  is: 


6  =  tan-1 


<  XA  ~XB  J 


(3.4) 


The  particle  and  their  angular  positions  for  tests  A  and  B  are  plotted  in  Figs.  35  and  36. 
Figure  26  shows  a  close-up  view  of  sample  A. 


Error  on  particle  angular  position 

As  shown  in  Fig.  34b,  the  error  on  particle  rotation  AO  can  be  calculated  from  the  error 
Ax  of  the  position  of  point  A  and  B. 


Ad  =  2sm'l(lAx/AB)^4Ax/AB  (3.5) 

where  is  it  assumed  that  Ax  «  AB.  The  error  AO  is  twice  as  large  for  the  smallest 
particles  (4  mm)  that  for  the  largest  particles  (8  mm). 


A' 

(a)  (b) 

Figure  34.  Definition  of  particle  rotation:  (a)  value  of  rotation  increment  and  (b)  error  on 
determination. 
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Figure  35.  Geometry  of  sample  A  at  four  different  stages  of  axial  compression  A1-A4. 
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Figure  36. 


Geometry  of  sample  B  at  four  different  stages  of  axial  compression  B1 


B4. 
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Interpretation  of  Results 

Once  the  center,  radius  and  rotation  of  particles  has  been  determined,  various  results  on 
the  kinematics  of  granular  materials  can  be  investigated  in  a  similar  way  to  previous 
numerical  investigations,  but  this  time  based  on  real  experimental  measurements  instead 
of  numerical  simulations. 

Stereophotogrammetric  visualization 

Figure  37  shows  a  three-dimensional  representation  of  the  false  relief  visualized  in  the 
stereograph  for  assemblies  of  discrete  particles.  The  rotation  of  particles  are  perceived  as 
if  the  particles  where  slanted  instead  of  being  cut  perpendicularly.  The  abrupt  jumps  in 
displacement  are  perceived  as  a  cliff.  Cluster  of  particles  undergoing  uniform  translation 
are  shown  as  a  flat  plateau,  while  rigid  clusters  under  rotation  are  represented  as  inclined 
plateau. 


Figure  37.  Displacement  and  rotation  of  particles  visualized  by  stereophotogrammetry. 


Particle  displacement 


As  shown  in  Figs.  38  and  39,  the  displacement  vectors  of  particle  centers  in  tests  A  and  B 
were  plotted  in  two  different  ways.  They  were  plotted  (1)  at  the  particles  centers,  which 
are  irregularly  spaced,  and  (2)  on  a  regularly  spaced  grid  by  interpolating  the 
displacements  vectors  at  the  particle  centers.  The  uniform  grid  displays  more  clearly 
strain  localization  (i.e.,  shear  band).  As  shown  in  Fig.  38  for  test  A,  a  shear  band  is 
formed  at  stages  A2-A3  and  A3-A4;  it  is  inclined  at  about  46°  with  respect  to  the 
horizontal.  As  shown  in  Fig.  39  for  test  B,  a  shear  band  is  displayed  at  stages  B2-B3  and 
B3-B4;  it  is  inclined  at  42  °  with  respect  to  the  horizontal.  Both  inclinations  are  smaller 
that  those  predicted  by  the  Mohr-Coulomb  theory,  i.e.,  9  =K/4  +  (p/2  when  <j>=180  for 
test  A  and  <|>=T  7  °  for  test  B. 

Shear  strain 

Figures  40  and  41  display  the  shear  strain  next  to  the  incremental  displacement  vectors  at 
the  particle  centers.  Shear  strains  were  computed  using  local  displacement  gradient 
centered  at  particle  centers,  as  described  in  Bardet  and  Proubet  (1991).  They  are 
represented  by  square  symbols  centered  on  the  particles;  their  sizes  are  proportional  to 
shear  strain  amplitude.  As  shown  in  Figs.  40  and  41,  shear  strains  are  concentrated  in 
narrow  band  that  coincide  with  the  shear  bands  of  Figs.  38  and  39. 

Particle  rotation 

Figures  42  and  43  display  the  particle  rotation  next  to  incremental  displacement  vectors. 
Rotation  angles  are  represented  by  a  pie  slice  centered  at  the  particle  center.  The  particles 
rotations  are  clearly  concentrated  in  the  shear  bands  as  the  shear  strain.  It  is  clear  that 
particle  rotation  concentrate  inside  shear  band.  This  result  agrees  with  those  of 
numerical  simulation  performed  by  Bardet  (1994).  Particle  rotations  have  a  determinant 
influence  on  the  localized  failure  mechanisms  of  idealized  granular  media. 
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Figure  39.  Interpolated  displacements  on  uniform  grid  (top  row)  and  displacements  at 
particle  centers  (bottom  row)  between  four  successive  stages  of  axial 
compression  B1-B4  of  Sample  B. 
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A3-A4 


Figure  40.  Incremental  shear  strains  (top  row)  and  incremental  displacements  of  particle 
centers  (bottom  row)  between  four  successive  stages  of  axial  compression 
A1-A4  of  Sample  A. 
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B1-B2 


B2-B3 


B3-B4 


B1-B2  B2-B3  B3-B4 


Figure  4 1 .  Incremental  shear  strains  (top  row)  and  incremental  displacements  of  particle 
centers  (bottom  row)  between  four  successive  stages  of  axial  compression 
B1-B4  of  Sample  B. 
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A1-A2  A2-A3  A3-A4 


A1-A2  A2-A3  A3-A4 


Figure  42.  Incremental  displacements  of  particle  centers  (top  row)  and  particle  rotations 
(bottom  row)  between  four  successive  stages  of  axial  compression  A1-A4  of 
Sample  A. 
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B1-B2  B2-B3  B3-B4 


Figure  43.  Incremental  displacements  of  particle  centers  (top  row)  and  particle  rotations 
(bottom  row)  between  four  successive  stages  of  axial  compression  B1-B4  of 
Sample  B. 

Discussion 

The  experimental  technique  presented  above  measures  successfully  the  displacement  and 
rotation  of  idealized  granular  material,  and  supports  some  of  the  results  from  past 
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numerical  analyses.  However,  the  present  procedure  is  time  consuming  for  large  number 
of  paticles  as  it  requires  digitizing  three  points  per  particle.  More  work  is  required  to 
develop  an  automatic  vision  system  that  retains  the  advantages  of  stereophotogrammetry. 


CONCLUSION 


# 


We  have  investigated  shear  band  instability  in  granular  media,  and  some  of  the  higher- 
order  continuum  theories  proposed  to  account  for  instability.  We  have  explored  the 
micro-macro  mechanics  transition  relevant  to  material  instability,  and  devised  laboratory 
experiments  for  measuring  the  kinematics  of  particles.  Our  methodology  combined 
computational  granular  mechanics,  higher-order  continuum  mechanics,  and  laboratory 
experiments. 

Our  comprehensive  review  of  past  work  on  computational  granular  mechanics  reveals  the 
diversity  of  numerical  codes  in  granular  mechanics,  and  the  large  number  of  assumptions 
entering  these  codes. 

Our  theoretical  investigation  shows  that  the  average  stress  in  granular  media,  as  defined 
from  virtual  work,  may  be  asymmetric  in  the  absence  of  contact  moments.  We  specify  the 
circumstances  and  amplitude  of  stress  asymmetry,  and  calculate  the  corresponding  couple 
stress  and  first  stress  moment.  We  also  show  that  the  average  stress  is  always  symmetric 
when  it  is  alternately  defined  by  using  statics  and  no  contact  moment.  The  stress 
asymmetry,  which  results  from  external  moments,  has  an  amplitude  that  decreases  with 
the  volume  size.  The  asymmetric  stress,  couple  stress  and  first  stress  moment  are 
analytically  calculated  in  particular  examples  relevant  to  shear  band  instability. 

Our  experimental  methodology  is  based  on  the  application  of  stereophotogrammetry  for 
measuring  the  kinematics  of  assembly  of  cylindrical  rods,  simulating  granular  media.  We 
were  able  to  measure  the  kinematics  of  a  large  number  of  nearly  identical  particles  in  the 
laboratory.  Our  experimental  methodology  has  generated  new  experimental  data  sets  on 
granular  media  useful  to  re-examine  the  applicability  of  higher-order  continua  on  the 
response  of  granular  media,  and  are  therefore  instrumental  to  understand  material 
instability  in  granular  media. 
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APPENDIX 


Listing  of  stereophotgrammetric  data  processing  programs  used  in  Part  III. 


GET_AB  calculates  of  particle  center  and  radius  by  three 
different  methods,  and  calculates  particle  rotation 

INPUT  FILES:  C1C2.0RI,C2C3.0RI,C3C4.0RI 

OUTPUT  FILES:  BALL_A.OUT,BALL_B.OUT,BALL_C.OUT 
ABOUT 

DIAMETER_A.OUT 
PARAMETER(M=2000) 

REAL  FX1(2,3),FY1(2,3),FX2(2,3),FY2(2,3),FX3(2,3),FY3(2,3) 

REAL  AX1(2,M),AY1(2,M),BX1(2,M),BY1(2,M),CX1(2,M),CY1(2,M), 

&  AX2(2,M),AY2(2,M),BX2(2,M),BY2(2,M),CX2(2,M),CY2(2,M), 

&  AX3(2,M),AY3(2,M),BX3(2,M),BY3(2,M),CX3(2,M),CY3(2,M) 

REAL  X1(2,M),Y1(2,M),R1(2,M),X2(2,M),Y2(2,M),R2(2,M), 

&  X3(2,M),Y3(2,M),R3(2,M) 

REAL  D0(M),TA(M),T1(M),T2(M),T3(M),T4(M),D1(2,M),D2(2,M),D3(2,M) 
CHARACTER*80  NAME,DIR1,DIR2,T 

C  Work  directory 

WRITE(*,'(A,$)') '  Enter  Test  A  or  B: ' 

READ(*/(A)')  T 

DIR1='F:\WORK\AFOSR\Processing\TEST//TRIM(T)//'\RawV 

DIR2=’F:\WORK\AFOSR\Processing\TEST7/TRIM(T)/AProcessedV 

C  Read  data  from  files 

NAME=TRIM(DIR1)//’C1C2.0RI' 

OPEN(1,  FILE=NAME,STATUS='OLD') 

CALL  READ_DATA(FX1  ,FY1  ,AX1  ,AY1  ,BX1  ,BY1  ,CX1  ,CY1  ,N1 , 1 ) 
CLOSE(I) 

NAME=TRIM(DIR1  )//’C2C3.0RI’ 

OPEN(1,  FILE=NAME,STATUS='OLD’) 

CALL  READ_DATA(FX2,FY2,AX2,AY2,BX2,BY2,CX2,CY2,N2,2) 
CLOSE(I) 

NAME=TRIM(DIR1)//’C3C4.0RI’ 

OPEN(1,  FILE=NAME,STATUS='OLD') 

CALL  READ_DATA(FX3,FY3,AX3,AY3,BX3,BY3,CX3,CY3,N3,2) 
CLOSE(I) 

C  Check  data 

WRITE(6,*)  'Fixed  points  1' 

WRITE(6,*)  'D12  =',DIST(FX1(1,1),FY1(1,1),FX1(1,2),FY1(1,2)) 
WRITE(6,*)  'D23  =,,DIST(FX1(1,3),FY1(1,3),FX1(1,2),FY1(1,2)) 
WRITE(6,*)  'Fixed  points  2' 

WRITE(6,*)  'D12  -  ,DIST(FX2(1,1),FY2(1,1),FX2(1,2),FY2(1,2)) 
WRITE(6,*)  'D23  =',DIST(FX2(1,3),FY2(1,3),FX2(1,2),FY2(1,2)) 
WRITE(6,*)  'Fixed  points  3' 

WRITE(6,*)  'D12  =',  DIST(FX3(  1 , 1  ),FY3(1 , 1 ),  FX3(  1 ,2),FY3(1 ,2)) 
WRITE(6,*)  'D23  =',DIST(FX3(1,3),FY3(1,3),FX3(1,2),FY3(1,2)) 
WRITE(6,*)  'Number  of  particles', N1  ,N2,N3 


C  Calculate  particle  rotation 
DO  1=1, N1 

T1  (l)=ATAN2(BY1  (1 ,1)-AY1(1 ,1),BX1  (1 ,1)-AX1  (1,1)) 
T2(I)=ATAN2(BY1(2,I)-AY1(2,I),BX1(2,I)-AX1(2,I)) 
T3(I)=ATAN2(BY3(1,I)-AY3(1,I),BX3(1,I)-AX3(1,I)) 
T4(I)=ATAN2(BY3(2,I)-AY3(2,I),BX3(2,I)-AX3(2,I)) 

END  DO 

C  Calculate  particle  center  and  radius  by  Method  A 

CALL  PARCEN_A(X1  ,Y1  ,R1  ,AX1  ,AY1  ,BX1  ,BY1  ,CX1  ,CY1  ,N1 ) 
CALL  PARCEN_A(X2,Y2,R2,AX2,AY2,BX2,BY2,CX2,CY2,N2) 
CALLPARCEN_A(X3,Y3,R3,AX3,AY3,BX3,BY3,CX3,CY3,N3) 

C  Write  ouput  results 

NAME=TRIM(DIR2)//’BALL_A.OUT 
OPEN(4,FILE=NAME, STATUS- UNKNOWN') 

CALL  WRITE_BALL(X1  ,Y1  ,R1  ,X3,Y3,R3,T1  ,T2,T3,T4,N1 ) 
CLOSE(4) 

NAME=TRIM(DIR2)//’DIAMETER_A.OUT 
OPEN(4,FILE=NAME1STATUS='UNKNOWNl) 
WRITE(4,'(6(5X,A4))')  'D1  ',,D2',,D2,,'D3',’D3',,D4' 

DO  1=1, N1 

WRITE(4,,(6F12.5)’)2.‘R1(:,I),2.*R2(:,I),2.*R3(:,I) 

END  DO 
CLOSE(4) 

C  Calculate  particle  center  and  radius  by  Method  B 
DO  1=1  ,N1 

AB=DIST(AX1(1,I),AY1(1,I).BX1(1,I),BY1(1,I)) 

D0(I)=(-(X1(1,I)-AX1(1,I))*(BY1(1,I)-AY1(1,I))+ 

&  (Y1(1,I)-AY1(1,I))*(BX1(1,I)-AX1(1,I))  )/AB 

END  DO 

CALL  PARCEN_B(X1  ,Y1  ,R1  ,AX1  ,AY1  ,BX1  ,BY1  ,D0,N1 ) 

CALL  PARCEN_B(X2,Y2,R2,AX2,AY2,  BX2,  BY2,  D0.N2) 

CALL  PARCEN_B(X3,Y3,R3,AX3,AY3,BX3,BY3,D0,N3) 

C  Write  ouput  results 

NAME=TRIM(DIR2)//’BALL_B.OUT 
OPEN(4,FILE=NAME, STATUS-UNKNOWN') 

CALL  WRITE_BALL(X1  ,Y1  ,R1  ,X3,Y3,R3,T1  ,T2,T3,T4,N1 ) 
CLOSE(4) 

C  Calculate  particle  center  and  radius  by  Method  C 
DO  1=1, N1 

TA(I)=D0(I)/DIST(AX1(1,I),AY1(1,I),BX1(1,I),BY1(1,I)) 

END  DO 

CALL  PARCEN_C(X1  ,Y1  ,R1  ,AX1  ,AY1  ,BX1  ,BY1  ,TA,N1 ) 

CALL  PARCEN_C(X2,Y2,R2,AX2,AY2,BX2,BY2,TA,N2) 

CALL  PARCEN_C(X3,Y3,R3,AX3,AY3,BX3,BY3,TA,N3) 

C  Write  ouput  results 

NAME=TRIM(DIR2)//’BALL_C.OUT 
OPEN(4,FILE=NAME, STATUS-UNKNOWN') 

CALL  WRITE_BALL(X1  ,Y1  ,R1  ,X3,Y3,  R3,T1  ,T2,T3,T4,N  1 ) 
CLOSE(4) 

C  Calculate  distance  AB  and  diameter 
DO  1=1  ,N1 
DO  J=1,2 

D1(J,I)=DIST(AX1(J,I),AY1(J,I),BX1(J,I),BY1(J,I)) 

D2(J,I)=DIST(AX2(J,I),AY2(J,I),BX2(J,I),BY2(J,I)) 

D3(J,I)=DIST(AX3(J,I),AY3(J,I),BX3(J,I),BY3(J,I)) 
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END  DO 
END  DO 

NAME=TRIM(DIR2)//7\B.OUT 
OPEN(4IFILE=NAMEISTATUS-UNKNOWN') 
WRITE(4,'(6(5X,A4))')  'AB1^,AB2,;AB2,;AB3,;AB3^,AB4, 


DO  1=1  ,N  1 

WRITE(4/(6F12.5)')  D1(:fl)ID2(:ll),D3(:ll) 

END  DO 

CLOSE(4) 

END 


SUBROUTINE  PARCENJ^X.Y.R.AX.AY.BX.BY.CX.CY.N) 

INTEGER  N 

REAL  X(2,N)1Y(2tN)IR(2IN)lAX(2lN),AY(21N)1  BX(2iN)IBY(2iN)iCX(2lN)1CY(2iN) 
C  PARCEN_A  calculate  the  center  and  radius  of  particles  with  3  points 
C  Input 

C  N:  total  number  of  particles 

C  AX.AY.BX.BY.CX.CY:  coordinates  of  points  A,  B  and  C 
C  Output 

C  X,Y:  center  position 

C  R:  radius 

DO  l=1,N 
DO  J=1,2 
A=AX(J,I)-BX(J,I) 

B=AY(J,I)-BY(J,I) 

E= (AY  ( J ,  I  )**2-B  Y  ( J ,  I  )**2 + AX(  J ,  I  )**2-BX(  J ,  I  )**2)/2. 

C=AX(J,I)-CX(J,I) 

D=AY(J,I)-CY(J,I) 

F=(AY(JIl)**2-CY(J,l)**2+AX(J,l)**2-CX(J1l)**2)/2. 

X(J,I)=(E*D-B*F)/(A*D-B*C) 

Y(Jtl)=(A*F-E*C)/(A*D-B*C) 

R(J1I)=DIST(X(J1I),Y(J1I),AX(JII),AY(J,I)) 

END  DO 
END  DO 
END 


SUBROUTINE  PARCEN^X.Y.R.AX.AY.BX^Y.DO.N) 

INTEGER  N 

REAL  X(2,N)iY(2IN),R(2tN)1AX(2,N),AY(2iN)fBX(2,N)iBY(2IN)lD0(N) 
C  PARCENJ3  calculate  the  center  and  radius  of  particles 
C  with  constant  vector  DO 
C  Input 

C  N:  total  number  of  particles 

C  AX,AYfBX,BY  coordinates  of  points  A  and  B 
C  DO  constant  length  DO 

C  Output 

C  X,Y:  center  position 

C  R:  radius 

DO  l=1,N 
DO  J=1,2 
XX=BX(J,I)-AX(J,I) 

YY=BY(J,I)-AY(J,I) 

X(J,l)=(AX(J,l)+BX(J,l))/2.-YY*D0(l)/SQRT(XX**2+YY**2) 
Y(J,l)=(AY(J,l)+BY(J,l))/2.+XX*D0(l)/SQRT(XX**2+YY**2) 
R(J,l)=DIST(X(Jf  l),Y(J,l),  AX(J,I),AY(J,I)) 

END  DO 
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END  DO 
END 


C====================================================== 

SUBROUTINE  PARCEN_C(X,Y,R, AX, AY, BX, BY, TA,N) 

INTEGER  N 

REAL  X(2,N),Y(2,N),R(2,N),AX(2,N),AY(2,N),BX(2,N),BY(2,N),TA(N) 
C  PARCEN_C  calculate  the  center  and  radius  of  particles  with  method  C 
C  with  constant  angle  ALPHA 
C  Input 

C  N:  total  number  of  particles 

C  AX, AY, BX, BY  coordinates  of  points  A  and  B 

C  ALPHA  constant  angle 

C  Output 

C  X,Y:  center  position 

C  R:  radius 

DO  1=1, N 
DO  J=1,2 

X(J,l)=(AX(J,l)+BX(J,l))/2.-TA(l)*(BY(J,l)-AY(J,l)) 

Y(J,l)=(AY(J,l)+BY(J,l))/2.+TA(l)*(BX(J,l)-AX(J,l)) 

R(J,I)=DIST(X(J,I),Y(J,I),AX(J,I),AY(J,I)) 

END  DO 
END  DO 
END 


SUBROUTI NE  READ_DATA(FX,  FY, AX, AY,  BX,  BY,CX,  CY,N,K) 
REAL  FX(2,1),FY(2,1), 

&  AX(2,1),AY(2,1),BX(2,1),BY(2,1),CX(2,1XCY(2L12 _ 

REAL  CNX(2,4),CNY(2,4) 

C  read  rotation  D  and  translation  (DX,DY)  of  global  data 
READ(1,*)  D1,DX1,DY1 
IF(K==2)  THEN 
READ(1,*)  D2,DX2,DY2 
ELSE 
D2=D1 
DX2=DX1 
DY2=DY1 
ENDIF 

C  read  XMIN,XMAX,YMIN,YMAX  of  global  data 
READ(1,*)  XMIN.XMAX.YMIN.YMAX 
C  read  fixed  points 
DO  1=1,3 

CALL  READ_PT(FX(1,I),FY(1,I),DX1,DY1,D1,DX2,DY2,D2,K) 
END  DO 

C  read  coordinates  of  corners 
DO  1=1,4 

CALL  READ_PT(CNX(1,I),CNY(1,I),DX1,DY1,D1,DX2,DY2,D2,K) 
END  DO 

C  read  all  points  A,  B  and  C 
C  and  select  origin  at  third  fixed  point 
READ  (1,*)N 
DO  1=1  ,N 

CALL  READ_PT(AX(1,I),AY(1,I),DX1,DY1,D1,DX2,DY2,D2,K) 
CALL  READ_PT(BX(1,I),BY(1,I),DX1,DY1,D1,DX2,DY2,D2,K) 
CALL  READ_PT(CX(1,I),CY(1,I),DX1,DY1,D1,DX2,DY2,D2,K) 
AX(:,I)=AX(:,I)-FX(:,3) 

AY(:,I)=AY(:,I)-FY(:,3) 


-  Page  119  - 


Qo  Qo  Qo 


BX(:,I)=BX(:,I)-FX(:,3) 

BY(:,l)=BY(:fl)-FY(:,3) 

CX(:,I)=CX(:,I)-FX(:I3) 

CY(:fl)=CY(:tl)-FY(:,3) 

END  DO 

CLOSE(I) 

END 


SUBROUTINE  READ_PT(X,Y,DX1  ,DY1,D1IDX21DY2ID2,K) 
REAL  X(2),Y(2) 

PARAMETER(S=2. 1 ) 

READ(1  ,*)  XI  ,Y1  tX2,Y2 
IF(K==1)  THEN 
X(1)=X1 
Y(1)=Y1 
ELSE 

X(1  )=(X1  -DX1  )*COS(D1  HY1 -DY1  )*SIN(D1 ) 

Y(1  )=(Y1-DY1  )*COS(D1  )+(X1-DX1  )*SIN(D1 ) 

ENDIF 

X(2)=(X2-DX2)*COS(D2MY2-DY2)*SIN(D2) 

Y(2)=(Y2-DY2)*COS(D2)+(X2-DX2)*SIN(D2) 

C  scale  length  to  obtain  mm 

x=s*x 

Y=S*Y 

END 


REAL  FUNCTION  DIST(X1,Y1  ,X2,Y2) 

DIST=SQRT((X1  -X2)**2+(Y  1  -Y2)**2) 
END 


SUBROUTINE  WRITE^ALLCXI.YI.RI.XS.YS.RS.TI^.TSJ^N) 
REAL  X1(2)N),Y1{2,N),R1(2,N)1X3(21N),Y3(2,N),R3(2,N), 

&  T1(N)IT2(N),T3(N),T4(N) 

WRITE(4,*)  N 

WRITE(4,’(16(5X,A5))')  'XI'/YI'/Rr/Tr/XZ/YZ/RZ/TZ, 

&  'X3':y3\'R3';t3':x4';y4';r4';t4' 

DO  1=1  ,N 

WRITE(41,(16(1X,F9.4))')X1(1,I)1Y1(1,I),R1(1,I),T1(I), 
X1(2,I)(Y1(2,I)1R1(2,I),T2(I), 

X3(1  J),Y3(1  fl),R3(1 ,1),T3(I), 
X3(2,I)IY3(2,I)JR3(2,I),T4(I) 

END  DO 
END 
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o  o  o 


DISP_ROT  plots  displacement  vector  and  rotation  of  particles 
INPUT  FILE:  BALL_A.OUT 

PARAMETER  (M=2000) 

CHARACTER*256  LABEL 

REAL  XI  (M),Y1(M),R1  (M),T  1  (M),X2(M),Y2(M),R2(M),T2(M), 

&  X3(M),Y3(M),R3(M),T3(M),X4(M),Y4(M),R4(M),T4(M) 

CHARACTER*80  NAME,DIR2,T 

SCALE=5. 

C  Work  directory 

WRITE(V(A,$)')  ’  Enter  Test  A  or  B: ' 

READ(*,'(A)')  T 

DIR2- F:\WORK\AFOSR\Processing\TEST//TRIM(T)/AProcessed\' 

C  Read  center  position  and  radius 
NAME=TRIM(DIR2)//’BALL_A.OUT 
OPEN(1, FI  LE=NAME, STATUS-OLD') 

READ(1,*)N1 
READ(1,'(A)')  LABEL 
DO  1=1, N1 

READ(1,*)  X1(I),Y1(I),R1(I),T1(I),X2(I),Y2(I),R2(I),T2(I), 

&  X3(I),Y3(I),R3(I),T3(I),X4(I),Y4(I),R4(I),T4(I) 

END  DO 
CLOSE(I) 

C  Plotting  extrema 

XM I  N=MI  N  VAL(X1  ( 1 :  N 1 )) 

XMAX=MAXVAL(X1  ( 1 :  N 1 )) 

YMIN=M!NVAL(Y1(1:N1)) 

YMAX=MAXVAL(Y  1  (1 :  N 1 )) 

W=ABS(XMIN-XMAX) 

H=ABS(YMI  N-YMAX) 

XMIN=XMIN-0.05*W 

XMAX=XMAX+0.05*W 

YMIN=YMIN-0.05*H 

YMAX=YMAX+0.05*H 

C  Draw  displacement  vectors 

NAME=TRIM(DIR2)//’DISP_ROT.GIF/GIF' 

CALL  PGBEGIN(0, NAME, 1,1) 

CALL  PGSCR(0,  1.,  1.,  1.)  Iblack  becomes  white 
CALL  PGSCR(1 ,  0.,  0.,  0.)  Iwhite  becomes  black 
CALL  PGQVP(0,XP1  ,XP2,YP1  ,YP2) 

CALL  PGVPORT(0. 1,0.3, 0.5, 0.9) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 

DO  1=1, N1 

CALL  ARROW(X1  (l),Y  1  (l),X2(l),Y2(l)) 

END  DO 

CALL  PGVPORT(0.3,0.5,0.5,0.9) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 

DO  1=1, N1 

CALL  ARROW(X2(l),Y2(l),X3(l),Y3(l)) 

END  DO 

CALL  PGVPORT(0.5,0.7,0.5,0.9) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 

DO  1=1, N1 

CALL  ARROW(X3(l),Y3(l),X4(l),Y4(l)) 

END  DO 
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C  Draw  rotation 

CALL  PGVPORT(0.1 ,0.3,0. 1 ,0.5) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 

DO  1=1, N1 

CALL  PIE(X1(I),Y1(I),R1(I),T1(I),T2(I), SCALE) 
END  DO 

CALL  PGVPORT(0.3,0.5,0.1,0.5) 

CALL  PGWNAD(XM!N,XMAX,YM!N,YMAX) 

DO  l=1,N1 

CALL  PIE(X2(I),Y2(I),R2(I),T2(I),T3(I),SCALE) 
END  DO 

CALL  PGVPORT(0.5,0.7,0.1,0.5) 

CALL  PGWNAD(XMiN,XMAX,YMIN,YMAX) 

DO  1=1, N1 

CALL  PIE(X3(1),Y3(I),R3(I),T3(I),T4(I), SCALE) 
END  DO 
CALL  PGEND 
END 


C*ARROW 

C+ 

SUBROUTINE  ARROW(XO,YO,X1,Y1) 

C 

C  ARROW  plots  a  vector  from  (xO.yO)  to  (xl  ,y1 ) 
C— 

REAL  X(3),Y(3) 

PARAMETER  (ANG=0.44,TIP=0.2,ZMIN=0.1) 
Z=SQRT((X1-X0)**2+(Y1-Y0)**2) 

IF(Z<ZMIN)  RETURN 

X(1)=X0 

Y(1)=Y0 

X(2)=X1 

Y(2)=Y1 

CALL  PGLINE(2,X,Y) 
ALF=ATAN2(Y1-Y0,X1-X0) 

X(1  )=X(2)-Z*COS(ALF-ANG)*TIP 
Y(1)=Y(2)-Z*SIN(ALF-ANG)*TIP 
X(3)=X(2)-Z*COS(ALF+ANG)*TI  P 
Y(3)=Y(2)-Z*SIN(ALF+ANG)*TIP 
CALL  PGLINE(3,X,Y) 

END 


C*P1E 

C+ 

SUBROUTINE  PIE(X0,Y0,R,T0,T1,S) 

C 

C  PIE  plot  a  pie  of  radius  R  at  (X0,Y0)  starting  from  angle  TO  to  T1 
C  with  a  scale  S 
C- 

PARAMETER  (N=11) 

REAL  X(N),Y(N) 

PI=ATAN(1.)*4. 

IF(ABS(T1  -T0)<0.001)  RETURN 

X(1)=X0 

Y(1)=Y0 

DT=S*(T1-TO)/FLOAT(N-2) 


IF(S*ABS(T1-TO)>2.*PI)DT=2.*PI/FLOAT(N-2) 
DO  1=2, N-1 
T=TO+FLOAT(I-1)*DT 
X(I)=X(1  )+R*COS(T) 

Y(I)=Y(1)+R*SIN(D 
END  DO 
X(N)=X(1) 

Y(N)=Y(1) 

CALL  PGSFS(I) 

CALL  PGPOLY(N,X,Y) 

END 


o  o  o  o  o 


STRAIN  plots  strain  field  (Desrues'  technique) 

INPUT  FILE:  BALL_A.OUT 

USE  MSIMSL 
PARAMETER  (M=2000) 

CHARACTER*256  LABEL 

REAL  X1(M),Y1(M),R1(M),T1(M),X2(M),Y2(M),R2(M),T2(M), 

&  X3(M),Y3(M),R3(M),T3(M),X4(M),Y4(M),R4(M),T4(M) 

CHARACTER*80  NAME.DIR2.T 

SCALE=5. 

C  Work  directory 

WRITE(7(At$y) '  Enter  Test  A  or  B: ' 

READ(V(A)')  T 

DIR2='F:\WORK\AFOSR\Processing\TESr//TRIM(T)/AProcessedV 


C  Read  center  position  and  radius 
NAME=TRIM(DIR2)//’BALL_A.OUT 
OPEN(1  ,FILE=NAME, STATUS- OLD’) 

READ(1,*)  N1 
READ(1,'(A)')  LABEL 
DO  l=1,N1 

READ(1  ,*)  XI  (l)f  Y1  (l),R1  (l).T1  (l)vX2(l)1Y2(l).R2(l),T2(l)fl 
&  X3(l),Y3(l)lR3(l)IT3(l),X4(l)iY4(l),R4(l)IT4(l) 

END  DO 
CLOSE(I) 

C  Plotting  extrema 

XMIN=MINVAL(X1(1:N1)) 

XMAX=MAXVAL(X1  (1 :  N 1 )) 

YMIN=MINVAL(Y1(1:N1)) 

YMAX=MAXVAL(  Y 1  (1  :N1 )) 

W=ABS(XM1N~XMAX) 

H=ABS(YM1N-YMAX) 

XM I N =XM !  N  -0. 05*  W 
XMAX=XMAX+0.05*W 
YMIN=YMIN-0.05*H 
YMAX=YMAX+0.05*H 

C  Plot  shear  strain  at  centers  of  grid 
NAME=TRIM(DIR2)//,Strain.gif/GIF’ 

CALL  PGBEGINCO.NAME.I.I) 

CALL  PGSCR(0,  1,1,1.)  Iblack  becomes  white 
CALL  PGSCR(1,  0.,  0.,  0.)  Iwhite  becomes  black 
CALL  PGQVP(0,XP1  ,XP2,YP1,YP2) 

CALL  PGVPORT(0. 1,0.3, 0.5, 0.9) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 

CALL  STRAIN(X1  ,Y1  ,X2,Y2,N1  .XMIN.XMAX.YMIN.YMAX) 
CALL  PGVPORT(0. 3,0. 5,0. 5,0.9) 

CALL  PGWNAD(XM!N,XMAX,YMIN,YMAX) 

CALL  STRAIN(X2,Y2,X3IY3,N1, XMIN.XMAX.YMIN.YMAX) 
CALL  PGVPORT(0.5,0.7,0.5,0.9) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 

CALL  STRAIN(X3,Y3,X4,Y4,N1  .XMIN.XMAX.YMIN.YMAX) 


C  Plot  displacement  vectors 

CALL  PGVPORT(0. 1,0.3, 0.1, 0.5) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 
DO  1=1  ,N1 


CALL  ARROW(X1(l),Y1(l),X2(l),Y2(l)) 

_  END  DO 

•  CALL  PGVPORT(0.3, 0.5, 0.1, 0.5) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 
DO  1=1, N1 

CALL  ARROW(X2(l),Y2(l),X3(l),Y3(l)) 

END  DO 

CALL  PGVPORT(0.5,0.7,0.1,0.5) 

CALL  PGWNAD(XMIN,XMAX,YMIN,YMAX) 

•  DO  1=1, N1 

CALL  ARROW(X3(l),Y3(l),X4(l),Y4(l)) 

END  DO 
CALL  PGEND 
END 


C*STRAIN 

C+ 

SUBROUTINE  STRAIN(X1,Y1,X2,Y2,N1,XMIN,XMAX,YMIN,YMAX) 

REAL  XI  (N1  ),Y1  (N 1  ),X2(N  1  ),Y2(N  1  ),XMIN,XMAX,  YMIN,  YMAX 
C 

C  STRAIN  plots  a  square  symbol  to  represent  shear  strain 
C— 

PARAMETER  (NX=25,NY=50) 

REAL  XL(2,4),VL(2,4),XG(NX),YG(NY),EPS(3) 

REAL  DET,SHL(3,4),SHG(3,4),B(3,8) 

REAL.ALLOCATABLE::  XY(:,:),WX(:),WY(:),VX(:,:),VY(:,:),GAM(:,:) 

C  Define  local  shape  functions 
CALL  QUAD_SHL(SHL) 

C  Define  grid  points 

XXMAX=XMAX-0. 1  *(XMAX-XMIN) 

XXMIN=XMIN+0.1*(XMAX-XMIN) 

DO  1=1, NX 

XG(I)=XXMI  N+(XXMAX-XXMI  N)*(l-1  )/FLOAT(NX-1 ) 

END  DO 

YYMAX=YMAX-0. 1  ‘(YMAX-YMIN) 

YYMIN=YMIN+0.1*(YMAX-YMIN) 

DO  1=1, NY 

YG(I)=YYMIN+(YYMAX-YYMIN)*(I-1  )/FLOAT(NY-1 ) 

END  DO 

C  Calculate  displacement  at  grid  points 

IF(ALLOCATED(XY))  DEALLOCATE(XY,WX,WY,VX,VY,GAM) 
ALLOCATE(XY(2,N1),WX(N1),WY(N1),VX(NX,NY),VY(NX,NY)) 
ALLOCATE(GAM(NX-1 ,  NY-1 )) 

XY(1,:)=X1(:N1) 

XY(2,:)=Y1(:N1) 

WX=X2(:N  1  )-X1  (:  N 1 ) 

WY=Y2(:  N 1  )-Y1  (:  N 1 ) 

CALL  SURF(N1,XY,WX,NX,NY,XG,YG,VX,NX) 

CALL  SURF(N1,XY,WY,NX,NY,XG,YG,VY,NX) 

C  Calculate  strain  at  grid  points 
DO  1=1, NX-1 
DO  J=1,NY-1 
XL(1,1)=XG(I) 

XL(1,2)=XG(I+1) 
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XL(1,3)=XG(I+1) 

XL(1,4)=XG(I) 

XL(2,1)=YG(J) 

XL(2t2)=YG(J) 

XL(2,3)=YG(J+1) 

XL(2,4)=YG(J+1) 

CALL  QUADSHG(XL,SHL,SHG,DET) 

B=0. 

B(1,1::2)=SHG(1i:) 

B(2,2::2)=SHG(2,:) 

B(3,1::2)=SHG(2,:) 

B(3,2::2)=SHG(1,:) 

VL(1,1)=VX(I,J) 

VL(1,2)=VX(I+1,J) 

VL(1,3)=VX(I+1,J) 

VL(1,4)=VX(I,J) 

VL(211)=VY(I,J) 

VL(2,2)=VY(I,J) 

VL(2t3)=VY(l,J+1) 

VL(2,4)=VY(I,J+1) 

EPS=MATMUL(B(:l:),RESHAPE(VL,(/8/))) 

GAM(I,J)=SQRT((EPS(1)-EPS(2))**2+EPS(3)**2) 

END  DO 
END  DO 

C  Plot  square  symbol  at  center  of  grid 
SMAX=10. 

GMAX=MAXVAL(GAM) 

DO  1=1,  NX-1 
DO  J=1,NY-1 
XX=(XG(I)+XG(I+1  ))/2. 

YY=(YG(J)+YG(J+1))/2. 

CALL  PLOT_SQUARE(XX,YY,GAM(l,J),GMAX,SI\/IAX) 
END  DO 
END  DO 
END 


SUBROUTINE  QUADSHG(XL,SHL,SHG,DET) 
REAL  XL(2,4),DET,SHL(314),SHG(3,4) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c= 


QUADSHG  calculates  global  derivatives  of  shape  functions  and 
jacobian  determinant  for  a  four  node  quadrilateral  element 
XL(J, I)  global  coordinate 

DET  jacobian  determinant 

SH(1 ,1)  local  XI  derivative  of  shape  function 
SH(2,I)  local  ETA  derivative  of  shape  function 

SH(3,I)  local  shape  function 

SHG(1  ,l)x  derivative  of  shape  function 
SHG(2, 1 )y  derivative  of  shape  function 
SHG(3,l)shape  function 

I  local  node  number  or  global  coordinate  number 

J  global  coordinate  number 

REAL  XS(2,2),TEMP 
SHG=SHL 

XS=MATMUL(SHG(1 :2,:),TRANSPOSE(XL(1 :2,:))) 

DET =XS(1 , 1  )*XS(2,2)-XS(  1 ,2)*XS(2, 1 ) 

XS=XS/DET 

DET=4.*DET 


DO  1=1,4 

TEMP=XS(2,2)*SHG(1 ,1)-XS(1 ,2)*SHG(2,I) 

SHG(2,I)=-XS(2,1)*SHG(1,I)+XS(1,1)*SHG(2,I) 

SHG(1,I)=TEMP 

END  DO 

END 


C*QUAD_SHL 

C+ 

SUBROUTINE  QUAD_SHL(SHL) 

REAL  SHL(3,4) 

C  QUAD_SHL  calculates  shape  functions  and  local  derivatives  for  element  QUAD 
C  SHL(1 ,1)  local  XI  derivative  of  shape  function 
C  SHL(2,I)  local  ETA  derivative  of  shape  function 
C  SHL(3,I)  local  shape  function 

C  I  local  node  number 

C— 

REAL  RG(4),SG(4) 

DATA  RG/-0.5,  0.5.0.5.-0.5/ 

DATA  SG/-0.5,-0.5,0.5,  0.5/ 

DO  1=1,4 

SHL(1,I)=0.5*RG(I) 

SHL(2,I)=0.5*SG(I) 

SHL(3,I)=0.25 
END  DO 
END 


C*PLOT_SQUARE 

C+ 

SUBROUTINE  PLOT_SQUARE(X,Y,V,VMAXtSMAX) 


C  PLOT_SQUARE  plots  a  square  at  position  (X,Y)  proportional  to  V 
C 

C  Input 

C  X.Y  location  where  to  center  square 
C  V  value 

C  VMAX  maximimum  value  corresponding  to  maximum  size  SMAX 
C  SMAX  maximum  size  for  square 
C 

c — 

REAL  XR(5),YR(5) 

C  define  square 
D= V/VM  AX*S  MAX 
XR(1)=X-D/2. 

YR(1)=Y-D/2. 

XR(2)=XR(1)+D 

YR(2)=YR(1) 

XR(3)=XR(2) 

YR(3)=YR(2)+D 

XR(4)=XR(3)-D 

YR(4)=YR(3) 

XR(5)=XR(1) 

YR(5)=YR(1) 

C  full  for  V>0,  hollow  for  V<0 


o  o  o  o  o  o  o 


CALL  PGSCI(I) 

CALL  PGSFS(2) 

CALL  PGP0LY(5,XR,YR) 

RETURN 

END 


‘ARROW 

+ 

SUBROUTINE  ARROW(XO,YO,X1,Y1) 
ARROW  plots  a  vector  from  (xO.yO)  to  (xl.yl) 


REAL  X(3),Y(3) 

PARAMETER  (ANG=0.44,TIP=0.2,ZMIN=0. 1 ) 
Z=SQRT((X1  -X0)**2+(Y1  -Y0)**2) 

IF(Z<ZMIN)  RETURN 

X(1)=X0 

Y(1)=Y0 

X(2)=X1 

Y(2)=Y1 

CALL  PGLINE(2,X,Y) 

ALF=ATAN2(Y1  -Y0.X1  -X0) 

X(1  )=X(2)-Z*COS(ALF-ANG)*TIP 
Y(1  )=Y(2)~Z*SIN(ALF-ANG)*TIP 
X(3)=X(2)-Z*COS(ALF+ANG)*TI  P 
Y(3)=Y(2)-Z*SIN(ALF+ANG)*T1P 
CALL  PGLINE(3,X,Y) 

END 


ACKNOWLEDGMENT 


This  document  reports  on  the  research  sponsored  by  the  AFOSR  grant  (F49620-93-1- 
0295)  augmented  by  the  AASERT  extension  (F49620-95-1-0420).  This  work  has  been 
partially  supported  by  the  Air  Force  Office  of  Scientific  Research,  the  National  Science 
Foundation,  and  by  the  ALERT  program  of  the  European  Community.  The  AFOSR  grant 
(F49620-93-1-0295)  was  augmented  by  an  AASERT  grant  (F49620-95-1-0420)  to 
support  Julie  Young,  who  subsequently  to  this  award  was  granted  a  NSF  fellowship  in 
1998.  The  authors  are  thankful  to  I.  Vardoulakis  of  the  National  Technical  University  of 
Athens,  Greece,  and  J.  Desrues  and  D.  Caillerie  of  the  University  Joseph  Fourier  of 
Grenoble,  France,  for  valuable  comments.  The  authors  also  thanks  Dr.  Nakase  of  the 
EPRI,  Japan,  for  making  available  its  discrete  materials. 


TABLE  OF  CONTENTS 

Abstract . 3 

Acknowledgment . 4 

Introduction . 1 

Shear  band  instability . 1 

Research  objectives  and  report  organization . 5 

Part  I.  Granular  mechanics . 6 

Introduction . ® 

Examples  of  discrete  modeling  in  soil  mechanics . 7 

Example  1 . 7 

Example  2 . 9 

Applications  of  discrete  modeling  in  engineering  and  applied  sciences . 10 

Discrete  and  continuous  modeling . 13 

Limitations  of  discrete  modeling  in  soil  mechanics . 14 

Numerical  methods  for  discrete  modeling . 14 

Physical  modeling  of  granular  media . 17 

Geometry  of  grains . 17 

Particles . 18 

Contacts . 20 

Contact  models . 31 

Governing  equations  of  statics . 35 

Boundary  conditions . 36 

Part  II.  Transition  from  discrete  to  continuous  media . 39 

Background . 39 

Granular  medium . 39 

Definition . 39 

Virtual  work  in  granular  media . 44 

Continuum  for  granular  media . 46 

Definition  of  average  stresses  in  granular  media . 48 

Average  stress . 48 


Symmetry  of  average  stress . 

Average  micropolar  stress  and  first  moment  of  stress . 

Alternate  definition  of  average  stress . 

Examples . 

Example  1 :  Double  layer  interface . 

Example  2:  Multi-layer  interface . 

Discussion . 

Conclusion . 

Part  III.  Experimental  investigation . 

Background . 

Axial  compression  of  idealized  granular  media . 

Sample  composition  and  fabrication . 

Experimental  setup  for  axial  compression . 

Results . 

Stereophotogrammetry  Measurement . 

Principle  of  stereophotogrammetry  for  measurement  of  displacement 

Application  of  stereophotogrammetry  to  kinematics  of  particles . 

Determination  of  Assembly  Geometry . 

Assumption  of  rigid  particles . 

Particle  center  and  radius . 

Particle  angular  position . 

Interpretation  of  Results . 

Stereophotogrammetric  visualization . 

Particle  displacement . 

Shear  strain . 

Particle  rotation . 

Discussion . 

Conclusion . 

References . 

Appendix . 


.49 

.50 

.50 

.52 

.52 

.54 

.58 

.58 

.59 

.59 

.61 

.61 

.63 

,64 

.65 

.65 

,66 

,68 

,68 

,69 

,75 

,78 

..78 

,79 

..79 

..79 

..85 

..87 

..88 

116 


INTRODUCTION 


Among  engineering  materials,  granular  materials  experience  the  most  complicated  and 
diversified  types  of  instabilities,  owing  to  their  particulate  and  multiphase  structures. 
These  instabilities  include  surface  instability  (e.g.,  rockburst  and  exfoliation),  volume 
instability  (e.g.,  buckling  and  arching),  localized  instability  (e.g.,  shear  bands),  and  fluid- 
grain  interaction  instability  (e.g.,  liquefaction).  These  various  types  of  instability  have 
been  observed  in  the  field,  in  association  with  global  and/or  local  failures  of 
geomaterials.  They  have  been  reproduced  and  studied  in  the  laboratoiy,  and  analyzed  by 
using  continuum  mechanics.  However,  the  physical  origins  of  instabilities  in  geomaterials 
are  still  poorly  understood. 

The  present  research  is  focused  on  shear  band  instability  in  granular  media,  and 
especially  on  the  physical  origins  of  these  instabilities.  In  the  following  review  of  past 
work,  we  intend  to  show  that  past  continuum  theories  are  based  on  diverse  assumptions, 
and  rely  on  semi-intuitive  arguments  not  necessarily  founded  on  physical  observations. 

Shear  band  instability 

Most  constitutive  models  for  geomaterials  describe  their  mechanical  behavior  through 
stress-strain  relationships,  which  have  no  internal  characteristic  length  specific  to 
materials.  This  macro-description  postulates  that  the  responses  of  small  laboratory 
specimens  are  the  scaled  responses  of  larger  granular  masses  in  the  field.  In  the  absence 
of  instability,  macro-continuum  mechanics  has  been  successfully  used  to  calculate 
stresses  and  strains  in  solids.  It  has  also  been  successful  to  explain  some  aspects  of  the 
unstable  behavior  of  geomaterials,  such  as  surface  and  volume  instability  (Hill  and 
Hutchinson,  1974;  Vardoulakis,  1979, 1981,  and  1988;  Vardoulakis  and  Graf,  1985; 
Bardet,  1991)  and  the  inclination  of  shear  bands  (Rudnicki  and  Rice,  1975;  Rice,  1976; 
Vardoulakis,  1980;  Bardet,  1991).  However,  without  an  internal  length,  macro-  theories 
simply  fail  to  describe  the  thickness  of  shear  bands,  which  was  found  to  be  about  eight 
times  the  mean  grain  size  in  granular  materials  (e.g.,  Milhlhaus  and  Vardoulakis,  1987). 
The  thickness  of  shear  bands  is  an  internal  length  that  does  not  scale  with  the  specimen 
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size.  When  a  continuum  theory  ignores  this  internal  length,  the  numerical  solution  of 
boundary  value  problems  leads  to  strain-softening  response  and  localized  deformation 
that  strongly  depend  on  mesh  sizes  (e.g.,  de  Borst  and  Sluys,  1992). 

Four  alternate  types  of  continuum  theory  with  an  internal  length  have  been  proposed  to 
describe  shear  band  instability  in  geomaterials.  These  approaches  are:  1)  higher  order 
strain  gradient  theory,  2)  micropolar  theory,  3)  nonlocal  theory,  and  4)  viscous 
regularization. 

Higher  order  strain  gradient  approach 

Zbib  and  Aifantis  (1989)  introduced  higher  order  strain  gradients  into  continuum  models 
to  ascribe  an  internal  length  to  materials  and  to  capture  the  thickness  of  shears  bands. 

Zbib  and  Aifantis  (1988)  derived  analytical  solutions  for  the  post-bifurcation  structures  of 
shear  bands,  which  are  in  agreement  with  experimental  observations  on  metals. 
Vardoulakis  and  Aifantis  (1989)  investigated  the  effects  of  gradient  dependent  dilatancy 
on  shear  bands  and  liquefaction  instabilities  within  saturated  sands.  Their  assumption  of 
gradient  dependent  dilatancy  invokes  basic  microscopic  considerations  for  two- 
dimensional  materials.  Their  theory  produces  the  thickness  of  shear  bands,  and  creates 
liquefying  strips  in  saturated  sand  specimens.  However,  in  their  dilatancy  relationship, 
they  replaced  a  factor  1/2  by  an  "unspecified  dimensionless  coefficient  to  be  determined 
from  experiment."  They  did  not  clearly  explain  the  physical  origins  of  gradient  dependent 
dilatancy  and  liquefaction  patterns  either. 

The  gradient  effects  on  yield  and  plastic  potential  surfaces  are  difficult  to  exhibit  and 
calibrate  by  experiments  on  geomaterials.  Oka  et  al.  (1991)  applied  the  higher  order 
gradient  approach  to  clays.  Oka  calibrated  the  coefficients  of  higher  order  effects  by  trial 
and  error  to  produce  the  desirable  shear  band  thickness.  There  is  a  definite  need  for 
developing  a  approach  for  investigating  the  physical  justification  of  higher  order 
gradients  in  particulate  media. 
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Micropolar  approach 

The  micropolar  theory  (Eringen,  1966, 1 968)  is  a  continuum  version  of  the  structural 
theory  of  Cosserat  (1909).  It  enriches  the  kinematics  and  kinetics  of  continua  by  adding 
material  rotations  and  couple  stresses.  Using  a  micropolar  approach  similar  to  Kanatani 
(1979),  Muhlhaus  and  Vardoulakis  (1986)  explained  the  emergence,  orientation  and 
thickness  of  shear  bands  in  granular  materials.  Bardet  and  Proubet  (1992a)  used  a  similar 
linear  stability  analysis  and  micropolar  description,  and  investigated  the  structure  of 
persistent  shear  bands  in  idealized  granular  media.  They  successfully  described  the 
thickness  of  shear  bands  and  the  relation  between  particle  rotation  and  displacement 
within  persistent  shear  bands.  However,  the  coefficients  of  their  micropolar  models, 
based  on  the  flow  or  deformation  theory  of  plasticity,  had  to  be  set  to  unrealistic  values  to 
reproduce  the  observations. 

Chang  et  al.  (1990,  1991, 1992)  developed  micropolar  theories  for  granular  materials 
based  on  microscopic  models.  They  derived  the  micropolar  constants  in  terms  of  the 
inter-particle  stiffness,  and  investigated  the  micropolar  effects  on  the  solution  of  selected 
boundary  value  problems.  Chang  derived  stress-strain  relationships  without  examining 
their  effects  on  material  instability.  He  did  not  investigate  the  problem  of  strain 
localization  as  De  Borst  and  Sluys  (1992). 

Dietsche  et  al  (1991)  examined  the  micropolar  effects  on  bifurcation,  without  justifying 
the  physical  origins  of  their  continuum  assumptions.  Their  study  assesses  only  the 
micropolar  effects  on  the  uniqueness  of  boundary  value  problems. 

Like  the  higher  order  gradient  theory,  the  micropolar  theory  introduces  an  internal  length 
in  boundary  value  problems.  However,  there  is  still  need  for  establishing  the  physical 
relevance  of  the  micropolar  theory  for  the  instability  in  particulate  media. 

Nonlocal  continuum  approach 

Nonlocal  continuum  mechanics  (Bazant,  1991  and  1992;  Eringen,  1992;  Valanis,  1992)  is 
especially  suited  for  investigating  material  instability  and  strain  softening.  Nonlocal 
continuum  mechanics  replaces  the  local  stress  by  an  average  stress,  which  is  the  integral 


of  the  weighted  local  stress  over  a  material  volume.  Nonlocal  models  possess  an  internal 
length  that  relates  to  the  thickness  of  shear  bands  (Bazant  and  Pijaudier-Cabot,  1988). 
Nonlocal  models  have  been  especially  designed  to  assess  damage  in  brittle  materials 
(e.g.,  concrete)  with  cohesive  granular  bonds  (Ju,  1991).  Recently,  Bazant  (1992) 
introduced  a  new  concept  of  nonlocal  damage  based  on  the  micromechanics  of  crack 
interaction,  and  therefore  established  the  physics  for  the  spatial  averaging  integral  and 
weight  function  employed  in  nonlocal  models.  However,  the  concept  of  crack  interaction 
and  propagation  does  not  apply  to  a  random  packing  of  particulate  media  with 
cohesionless  contacts.  The  damage  concept  seems  more  applicable  to  cemented  or 
overconsolidated  soils,  for  which  cohesive  bonds  can  be  damaged  and  broken.  However, 
the  damage  concept  appeals  for  describing  the  strain  softening  of  dense  sands,  which  is 
attributed  to  stress-dilatancy. 

Viscous  regularization 

The  dissipative  properties  of  physical  viscosity  introduce  an  internal  length,  which  relates 
to  the  thickness  of  shear  bands  in  metals  under  dynamic  loads  (e.g.,  Perzyna,  1991). 
Artificial  viscosity,  referred  to  as  viscous  regularization,  was  first  introduced  for 
capturing  shock  waves  (Lapidus,  1967).  Loret  and  Prevost  (1990)  adapted  it  for 
geomaterials  to  capture  the  thickness  of  shear  bands.  However,  the  viscosity  measured  for 
granular  materials  in  the  laboratory  is  much  smaller  than  the  artificial  viscosity  required 
for  capturing  the  thickness  of  shear  bands.  The  diffusion  of  pore-water  generates  viscous 
effects  (Bardet,  1992)  that  are  also  too  small  for  explaining  shear  bands.  The  introduction 
of  artificial  viscosity  in  geomaterials  seems  to  be  a  numerical  expedient  deprived  of 
physical  justification.  This  numerical  artifice  requires  an  excessively  large  artificial 
viscosity  which  unrealistically  attenuates  waves  and  accelerations  within  boundary  value 
problems. 

The  continuum  theories  -  higher  order  gradient  theory,  micropolar  theory,  nonlocal 
theory,  and  viscous  regularization  -  agree  on  the  necessity  of  introducing  an  internal 
length.  They  justify  their  assumptions  by  the  need  to  limit  localization  in  the  numerical 
solution  of  boundary  value  problems.  Their  diverse  assumptions  rely  on  semi-intuitive 
arguments  not  necessarily  founded  on  microscopic  observations  for  particulate  media. 
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